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ABSTRACT 


The  use  of  the  digital  computer  for  the  simultaneous  generation  of 
circuit  responses  and  their  corresponding  sensitivity  to  parameter  changes 
is  considered  for  circuits  containing  diodes  and  transistors.  Modeling  of 
the  diodes  and  transistors  is  based  upon  the  Ebers-Moll  equations,  with 
the  addition  of  voltage-dependent  capacitors  to  account  for  switching 
timeo  The  resulting  general  iterative  equations  are  then  used  to  calcu- 
late the  effects  of  radiation  on  a  p-n  junction.  The  waveforms  and 
circuit  recovery  time  are  studied  as  a  function  of  circuit  parameters. 
The  use  of  sensitivity  functions  for  predicting  waveform  variations  is 
considered. 
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I.   INTRODUCTION 

A  measure  of  the  sensitivity  of  circuit-response  variables  to  param- 
eter variations  is  useful  in  circuit  analysis  and  design.  Many  different 
definitions  of  sensitivity  are  available.  H.  Bode  [1]  defines  a  logarith- 
mic or  normalized  sensitivity  function  of  a  variable,  x.  (t,a.),  with 

x.   3  In  a.  a 

respect  to  a  parameter,  a . ,  as  S   =   ,   J  .  The  inverse  of  this 

definition  is  adopted  by  I.  M.  Horowitz  [2].  Calculation  of  the  sensi- 
tivity of  linear-circuit  responses  to  incremental  changes  in  network 
parameters  has  been  presented  by  J.  Leeds  [3].  Recently  S.  R.  Parker  [4] 

generalized  the  technique  to  include  nonlinear  circuits.  Both  Leeds  and 

x.   3X. 
Parker  define  the  sensitivity  function  as  Sa  =  r— -  .  This  work  was 

aj   3aj 
preceded  by  that  of  P.  Kokotovic  [5]  and  R.  Tomovic  [6]  who  considered 

variations  in  the  solution  of  a  linear  differential  equation  due  to  incre- 
mental changes  in  the  coefficients.  Tomovic  defines  the  sensitivity 
x.   3x.  ,  x. 

coefficients  as  S  ]  =  -r—   ,  and  Kokotovic  adopts  the  definition  S  n  = 
sy  a .   aa .  a. 

3*,-  J      J  J 

1     ,  and  discusses  implementation  of  the  sensitivity  coefficient 
using  analog  computation.  The  work  by  Leeds  and  Parker  is  oriented 
towards  the  digital  computer.  In  this  thesis  digital  calculation  of 
sensitivity  functions  of  nonlinear  circuits  is  considered  in  detail  with 
specific  application  to  the  study  of  the  effects  of  radiation  on  solid- 

X.    3X. 

state  junctions.  The  sensitivity  function  is  defined  as  S   =  — ,     = 

ax.  aj   3  ln  aj 

/   ,  k    .  This  definition  is  adopted  since  it  enables  sensitivity  func- 

j  J 
tions  with  respect  to  parameters  having  widely  different  nominal  values  to 

be  conveniently  compared  on  a  percentage  basis,  and  also  avoids  some 

difficulties  encountered  in  the  definitions  adopted  by  Bode  and  Horowitz 

when  the  circuit-response  variable  is  zero.  It  should  be  noted  that  the 
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above  definitions  of  the  sensitivity  function  apply  only  for  infinites- 
imally  small  changes  in  x.  and  a-..  The  practicality  of  this  definition 
when  applied  to  specific  examples  with  large  parameter  increments  is 
considered  in  this  thesis. 

In  Chapter  II  it  is  shown  that  for  linear  circuits  the  sensitivity 
function  may  be  derived  using  the  compensation  theorem  [7]„  A  convenient 
set  of  iterative  equations  is  developed  using  the  concept  of  state  vari- 
ables which  allows  the  time  response  of  the  states  and  their  corresponding 
sensitivities  to  be  obtained  simultaneous ly„  A  linear  example  illustrates 
the  procedure . 

In  order  to  obtain  useful  results  from  the  sophisticated  circuit  anal- 
ysis and  design  programs  for  nonlinear  circuits,  acceptable  nonlinear 
models  are  required„  Chapter  III  discusses  modeling  of  diodes  and  tran- 
sistors for  digital  computer  analysis  based  on  the  Ebers-Moll  [8]  equations 
These  models  contain  nonlinear  capacitances  and  nonlinear  resistances  for 
the  p-n  junction   A  state-variable  formulation  is  developed  which  results 
in  a  set  of  iterative  equations  whose  solution  yields  time  responses  for 
circuits  containing  linear  elements,  diodes,  and  transistors .  The  proce- 
dure is  illustrated  by  considering  the  effect  of  a  radiation  current  pulse 
injected  across  (1)  a  p-n  junction  diode,  and  (2)  the  base-emitter  p-n 
junction  of  a  transistor  operating  in  the  active  region.  These  examples 
yield  new  results  concerning  the  effect  of  radiation  on  p-n  junctions  and 
comparative  data  for  several  runs  is  presented „ 

Sensitivity  analysis  of  nonlinear  circuits  is  discussed  in  Chapter  IV. 
The  state- variable  formulation  developed  in  Chapter  III  is  extended,  re- 
sulting in  a  set  of  iterative  equations  which  allow  time  responses  and 
their  corresponding  sensitivities  to  be  obtained  simultaneously.  These 
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procedures  are  applied  to  the  diode  example  of  the  previous  chapter  and 
the  results  interpreted  in  terms  of  the  sensitivity  of  the  diode  response 
to  percentage  parameter  variations. 
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II.  SENSITIVITY  ANALYSIS  OF  LINEAR  CIRCUITS 

A.  INTRODUCTION 

There  are  several  general  theorems  which  are  useful  in  analysis  and 
design  of  linear  networks.  One  of  these,  the  compensation  theorem  [7], 
has  been  known  for  many  years  and  may  be  applied  when  it  is  desired  to 
determine  what  effect  an  incremental  impedance  change  in  one  branch  of  a 
network  has  upon  the  currents  and  voltages  of  other  branches  of  the  net- 
work. Recently  the  concept  of  the  sensitivity  of  certain  circuit-response 
variables  to  incremental  changes  in  given  network  parameters  has  been 
presented  by  J,  Leeds  [3].  This  work  was  preceded  by  that  of  P.  Kokotovic 
[5]  and  R.  Tomovic  [6]  who  considered  variations  in  the  solution  of  a 
linear  differential  equation  due  to  incremental  changes  in  the  coeffi-. 
cientSc  Tomovic  and  Kokotovic  implemented  the  sensitivity  coefficient 
by  using  analog  computation,  Leeds'  work  is  oriented  towards  the  digital 
computer. 

In  the  next  section  sensitivity  analysis  of  linear  circuits  is  con- 
sidered. In  the  succeeding  section  the  relationship  between  the  older 
compensation  theorem  and  the  recent  sensitivity  concept  is  discussed. 
This  correlation  has  not  appeared  in  the  literature  heretofore.  The  final 
section  contains  an  illustrative  example  to  demonstrate  the  implementation 
of  the  results  on  a  digital  computer. 

B.  DETERMINATION  OF  SENSITIVITY  FUNCTIONS 
1 .  Auxiliary  Network  Approach 

Leeds  [3]  developed  a  relationship  between  the  sensitivity  func- 

xi    3xi 
tion,  S   =       (where  x.  is  a  response  variable  and  a.  is  a  circuit 

aj    3aj         "  J 
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parameter),  and  node  voltages  and  branch  currents  of  an  auxiliary  network 
called  the  sensitivity  model.  The  auxiliary  network  is  obtained  by  re- 
ducing all  independent  sources  to  zero  and  adding  a  single  source  in  the 
branch  containing  the  element  under  consideration.  The  source  depends 
upon  the  nature  of  the  element.  Fig.  2.1  summarizes  the  sources  to  be 
added  for  branches  containing  resistance,  inductance,  and  capacitance 
respectively.  The  voltages  and  currents  in  the  auxiliary  network  are  the 
voltage  and  current  sensitivities  with  respect  to  the  parameter  being 
considered. 

The  advantage  of  Leeds'  approach  is  that  it  allows  the  sensitivity  of 
a  network  to  be  calculated  by  evaluating  the  responses  of  a  coupled  auxil- 
iary network  rather  than  actually  performing  the  differentiation  indicated 
by  the  definition  of  the  sensitivity  function. 

2.  Augmented  State  Equations 

Sensitivity  functions  may  also  be  obtained  using  the  concepts  of 
state  variables.     The  network  is  first  represented  in  state-variable  form 
as  discussed  by  T.   Bashkow  [9],  and  Kuh  and  Rohrer  [10]. 

x     =     Ax  +  Bu_  (2.1) 

where 

x_  is  an  n  vector  of  circuit  states 
A  is  a  constant  nxn  matrix 

u_  is  an  m  vector  of  forcing  functions. 

x.     ax. 

The  sensitivity  functions,  S   =  r-5 ,  may  be  obtained  by 

a  •   o  n  a  . 
J        J 
performing  the  indicated  partial  differentiation  on  (2.1).  The  resulting 

partial  derivatives  can  be  considered  as  additional  states  which  may  be 

solved  in  conjunction  with  the  original  state  equations.  They  may  be 

included  in  the  following  augmented  matrix; 
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ORIGINAL  CIRCUIT  BRANCH  AUXILIARY   CIRCUIT   BRANCH 
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FIG.  2.1.     CIRCUIT    PARAMETER    BRANCHES    AND 
SENSITIVITY    MODEL   BRANCHES. 
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where     j  =  l,2,...,r   ,  and  r  is  the  number  of  sensitivity  parameters  and 
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If  u_  is  independent  of  a.,  that  is  if  a.  is  a  circuit  parameter 

J  J 


only,  (2.3a)  takes  the  form 
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Separating  the  set  of  first-order  partial   differential   equations 
of  (2,3a)   for  each  sensitivity  parameter  yields 

(2.5) 


-sj 
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+     A    .  x     +     B    .   u 
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for  j  =  1,2, 


Trapezoidal  integration  is  commonly  used  to  solve  a  set  of  first- 
order  matrix  differential  equations  because  the  procedure  always  converges 
and  provides  acceptable  accuracy  for  most  problems  [11].  Utilizing  trape- 
zoidal integration,  (2.1)  and  (2.5)  yield  the  following  iterative  solutions 
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where 


4>(T)     ■ 
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Equations   (2.6)  and  (2.7)  are  a  set  of  iterative  equations  whose 
solution  yields  the  time  response  of  the  states  and  their  corresponding 
sensitivities.     Combining  the  results  into  a  single  equation  results  in  a 
solution  format  similar  to  (2.3). 
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C.   USE  OF  THE  COMPENSATION  THEOREM  TO  DERIVE  SENSITIVITY  FUNCTIONS 

The  compensation  theorem  relates  an  incremental  impedance  change  in 
one  branch  of  a  network  with  corresponding  incremental  current  and  voltage 
changes  in  the  other  branches  [7],  In  order  to  relate  the  compensation 
theorem  directly  with  sensitivity  analysis  an  analagous  relation  for  incre- 
mental changes  in  the  branch  parameters  R,  L,  and  C  is  developed. 

1 .  Resistive  Change 

If  the  resistance  in  a  branch  is  changed  by  an  amount  AR,  and  a 
voltage  source  equal  to  I°(aR)  is  introduced  (Fig,  2,2a)  the  network  vari- 
ables remain  unchanged.  If  another  source  is  introduced  (Fig.  2.2b)  the 
branch  voltage  and  current  will  change  by  an  incremental  amount  aV  and  Al 
respectively.  By  applying  the  superposition  theorem  (Fig.  2.2c)  the 
incremental  voltage  and  current  can  be  calculated, 

AV  =  Al'(R  +  ARj  +  I -(AR)  (2.9) 

Dividing  (2,9)  by  aR  and  letting  ar^O  results  in  (2.10), 

«.««♦!  (2.10) 

2.  Inductive  Change 

If  the  inductance  in  a  branch  is  changed  by  an  amount  aL,  and  a 

A  T 

voltage  source  equal  to  aL  *  {-rr)   is  introduced  (Fig,  2,3a)  the  network 
variables  remain  unchanged.  If  another  source  is  introduced  (Fig.  2.3b) 
the  branch  voltage  and  current  will  change  by  an  incremental  amount  aV  and 
aI  respectively.  By  applying  the  superposition  theorem  (Fig.  2,3c)  the 
incremental  current  and  voltage  can  be  calculated. 


AV  =  (L  +  AL)  d|l  +  AL^i  (2.11) 
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Dividing  by  aL  and  letting  aL+O  results  in  (2,12), 

3L  L  dt  lal;       dt  ^-'^ 

3.  Capacitive  Change 

If  the  capacitance  in  a  branch  is  changed  by  an  amount  aC,  and  a 
current  source  equal  to  aO  (-nr)  is  introduced  (Fig.  2,4a)  the  network 
variables  remain  unchanged.  Introducing  another  current  source  (Fig.  2.4b) 
results  in  incremental  changes  in  the  branch  voltage  and  current.  Again 
by  applying  the  superposition  theorem  (Fig.  2.4c)  the  incremental  voltage 
and  current  can  be  calculated, 

Al  =  (C+  AC)  %f  +   aC^-  (2.13) 

Dividing  (2.13)  by  aC  and  letting  aC+0  results  in  (2.14). 

3C     L  dt  [dC}       dt  Uel4; 

Equations  (2.10),  (2.12),  and  (2.14)  are  identical  to  the  auxil- 
iary network  equations  of  Fig.  2.1.  This  illustrates  the  direct 
relationship  between  the  compensation  theorem  and  sensitivity  analysis. 

D.   LINEAR  EXAMPLE 

Consider  the  peaking  circuit  [12]  in  Fig,  2,5a  operating  in  the  linear 
region.  The  incremental  equivalent  circuit  is  represented  in  Fig.  2.5b 
where  R,  =  r  is  the  plate  resistance  and  C  is  the  coil  capacitance  plus 

output  capacitance,  stray  capacitance  and  wiring  capacitance.  The  equiv- 

R1R2 
alent  circuit  is  indicated  in  Fig.  2.5c  where  R  =  D-  -  -  D   and  v  = 

R^2  Kl   K2 

"  R]  +  R2  yVi  * 
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It  is  of  interest  to  determine  how  sensitive  the  output  is  with 
respect  to  incremental  changes  in  R,  L  and  C   As  representative  values 
let  L  =  0„1  henry,  C  =  250yf arads ,  and  R  be  adjustable  according  to  the 
type  of  response  desired.  The  response  can  be  determined  from  the  differ- 
ential equations  of  the  circuit  and  is  underdamped  for  K <  1  (R  =  20), 
critically  damped  for  K  =  1  (R  -   10),  and  overdamped  for  K  >  1  (R  =  5); 
where  K  =  or/L/C  . 

The  sensitivity  functions  may  be  obtained  by  constructing  auxiliary 
networks  as  outlined  in  the  preceding  section,  or  by  differentiating  (2.1) 
with  respect  to  In  R,  In  L,  and  In  C  respectively. 

The  state-variable  representation  is  illustrated  below.  The  solution 
takes  the  form  of  (2.8)  where 

u  =  v  =  input  voltage  of  70  volts  (2.15a) 


x  = 


output  voltage 
inductor  current 


A  - 


-1/RC 
1/L 


-1/C 


o  J  • 


B  - 


1/RC 
0 


(2.15b) 


(2.15c) 


usj  =  0   for  j  -  1,2,3 


(2.15d) 


-si 


3X, 

3  In  R 

3X2 

_3  In  R 

sensitivity  of  output  with  respect 
to  R 


sensitivity  of  inductor  current  with 
respect  to  R 


(2.15e) 


-S2 


3X, 

| 

3  In  L 

3Xp 

3  In  L 

sensitivity  of  output  with  respect 
to  L 

sensitivity  of  inductor  current  with 
respect  to  L 


(2.15f) 
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'si 


■1/RC 
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'S3 
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-1/RC 
0 


(2.15i) 


The  computer  program  used  to  implement  the  solution  is  included  in 
Appendix  B.     This  program  is  general   in  that  it  is  easily  modified  to 
handle  linear  circuits  containing  more  state  variables  and  sensitivity 
parameters. 

The  time  response  of  the  output  of  the  circuit  and  the  corresponding 
sensitivities  are  presented  in  Figs.   2.6  through  2.11   for  underdamped, 
critically  damped  and  overdamped  responses.     These  results  have  been 
obtained  using  trapezoidal   integration  as  well   as  a  fourth-order  Runge- 
Kutta  integration  technique. 
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Ill 0  MODELING  OF  NONLINEAR  CIRCUITS 

A.  INTRODUCTION 

The  digital  computer  has  played  an  increasingly  important  role  in 
circuit  analysis  and  design  in  the  past  ten  years ,  This  has  been  partly 
due  to  the  emergence  of  sophisticated  circuit-analysis  and  design  programs 
[13].  In  order  for  these  programs  to  be  useful,  acceptable  nonlinear 
models  which  relate  model  parameters  to  physical  processes  are  required. 
These  models  must  lend  themselves  to  numerical  analysis  and  provide  accept- 
able quantitative  accuracy „ 

In  the  next  section,  a  general  mathematical  format  which  is  convenient 
for  sensitivity  analysis  of  nonlinear  circuits  is  developed.  In  the  suc- 
ceeding section  modeling  of  diodes  and  transistors  is  outlined.  This  is 
followed  by  examples  which  consider  calculation  of  the  response  of  a 
p-n  junction  to  a  pulse  of  radiation  current, 

B,  GENERAL  NONLINEAR  STATE-VARIABLE  FORMAT 

The  state  equation  for  a  nonlinear  time-invariant  system  is  of  the  form 

x(t)     =     a(x(t),u_(t))  (3.1) 

where  a  is  a  vector  of  nonlinear  functions  of  the  state  variables  and 
forcing  functions. 

As  is  shown  in  the  section  which  follows,  (3.1)  takes  the  following 
general  form  for  circuits  containing  diodes  and  transistors  (represented 
by  the  Ebers-Moll  [8,14]  equations).  This  representation  includes  non- 
linear resistances  and  nonlinear  capacitances. 

x(t)  =  A(x(t)).x(t)  +  c(x(t))  +  B(x(t))-u(t)  (3.2) 
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where  A(x_(t))  is  a  nxn  matrix  and  B(_x(t))  is  an  nxm  matrix,  both  with 
elements  which  may  be  nonlinear  functions  of  the  states,  and  c_(x_(t))  is 
a  general  nonlinear  vector. 

Equation  (3.2)  may  be  written  as 

k{t)     =     f(x(t))  +  B(x(t)).u(t)  (3.3) 

where 

f(x(t))  =  A(x(t))-x(t)  +  c(x(t)) 

Using  trapezoidal   integration  for  (3.3)   results  in 

x(n)-x(n-l)     =     4jr  [f(x(n))  +  f(x(n-l))] 

+     £|  [B(x(n)).u(n)     +     B(x(n-1 ) ).  u_(n-l )]  (3.4) 

Equation  (3.4)  is  a  nonlinear  matrix  equation  implicit  in  xjn)  which, 
as  a  result,  is  generally  difficult  to  solve.  Simplification  results  if 
f_(x_(n))  and  B(x_(n) )•  u_(n)  ,  where  u_(n)  is  considered  constant,  are  expanded 
in  a  Taylor  series  about  the  point  x(n-l).  In  matrix  notation 


f(x(n))  =   z  (AxTv)k  f(x(n-l))/k!  (3.5a) 

k=0 

B(x(n)).u(n)  =  I     (AxTv)k  [B(x(n-1 )). u(n)]/k!  (3.5b) 

k=0 

where 

v  =  [  -i-  -i-  ...   4-  ]T 

3Xn    oX0  9X„ 

I    d  n 

f,  Bu_,  x_,  ax_  are  n  vectors 

and 

ax  =  x(n)  -  x(n-l ) 
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A  first-order  approximation  for  (3.5)   results  in  (3.6). 

f(x(n))     -     f(x(n-l))     +     (axTv)  I(x(n-1))  (3.6a) 

B(x(n)).u(n-)     =     B(x(n-1  ))«u(n)     +     (axTv)    [B(x(n-1 )).  u(n)]   (3.6b) 

Substituting   (3,6)   into  (3.4)  yields 

x(n)  -  x(n-l)     -    ^  [  2f(x(n-l))  +  (axTv)  f(x(n-l))   ] 

+    ^  [B(x(n-1))   •    (u(n)  +  u(n-l))  +  (axTv)   {B(x(n-1  ))-u(n)}] 

(3.7) 

Equation  (3.7)  may  be  simplified  by  using  the  matrix  identity 

(axTv)  f  =  (vfT)T  ax  (3.8) 

which  is  proven  in  Appendix  A. 

Substituting  (3.8)  into  (3.7)  and  solving  for  x_(n)  results  in 

x(n)     -     x(n-l)     +     aT- [I  -  ^-[(vfT)T  +  i  vOu(n)]T}T]  l"1 

[f  +  1b   ■    (u_(n)  +  u(n-l))]  (3.9) 

where  f,  B,  vf  and  (B-u_(n))   are  evaluated  at  the  point  xjn-1). 
Equation  (3.9)  can  be  solved  directly  for  x(n)  from  the  known  values  of 
x/n-1),  u_(n),  and  u^(n-1 ) . 

C.   DIODE  AND  TRANSISTOR  MODELING 

1 .  Modeling  for  a  p-n  Junction  Diode 

Diode  modeling  has  received  much  attention  during  the  past  twenty 
years  [15,16].  Recently  Steele  [17]  prepared  a  report  concerning  semi- 
conductor modeling  specifically  for  computer  analysis.  The  applicable 
results  are  summarized  in  this  section, 
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a.  Current-Voltage  Relationship 

The  current-voltage  relationship  for  an  ideal   diode  is 

id     =     rs  (eaV  -  1)  (3.10) 

where      I-  =  reverse  saturation  current 

a   =  a  parameter  dependent  upon  the  temperature  and  the 
type  of  semiconductor  material. 

Generally  the  reverse  saturation  current  is  not  entirely 
independent  of  voltage.  This  is  due  to  a  leakage  current  and  can  be 
accounted  for  by  the  addition  of  a  resistance,  R^,,,  shunting  the  diode. 
Also,  at  large  reverse  voltages  an  effect  known  as  Zener  breakdown  occurs.. 
This  effect  is  omitted  in  the  following  analysis. 

b.  Diode  Resistance 

The  static  resistance  Rj(v)  of  a  diode  is  defined  as  the  ratio 

v/ij.  Since  Rj-(v)  varies  widely  with  v  and  i  ,,  it  is  generally  not  a 

useful  parameter.  An  incremental,  or  dynamic,  resistance  r.(v)  is  defined 

as  r^(v)  =  -n—  = .  Although  the  dynamic  resistance  varies 

d  Eeote 

widely  with  v  it  is  a  useful   parameter  for  small-signal   operation. 

c.  Capacitance  Effect 

Diode  capacitance  is  due  to  two  separate  effects;  space  charge 
and  charge  storage. 

(1)  Space  Charge.  A  diffusion  potential,  Vn,  is  established 
in  the  junction  region  when  p  and  n  materials  are  placed  in  contact.  This: 
is  caused  by  an  increase  in  donor  atoms  on  the  n  side  and  acceptor  atoms 
on  the  p  side  of  the  junction,  resulting  in  a  capacitive  effect.  As  a 
reverse  voltage  is  applied  more  donor  and  acceptor  atoms  are  uncovered 
causing  the  junction  width  to  expand,  and  hence  the  capacitance  to  decrease, 
This  voltage-dependent  capacitance  has  the  following  form  when  the  diode 
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is  reverse  biased.  For  radiation  analysis,  acceptable  results  are  ob- 
tained when  Cn  is  assumed  to  be  constant, 

C0 
Cn  =  — =■  (3.11) 

D     (1  -  v/vD)n 

where     n  =  1/2  for  abrupt  junctions 

n  =  1/3  for  linear-graded  junctions. 
(2)  Charge  Storage,  When  the  diode  junction  is  forward  biased 
the  potential  barrier  is  lowered.  This  results  in  a  high  density  of  minor- 
ity carriers  injected  into  the  p  and  n  regions,  which  can  be  viewed  as 
stored  carriers.  A  change  in  voltage  results  in  a  change  in  density  (num- 
ber of  stored  carriers)  and  hence  a  capacitive  effect  . 

If  the  voltage  on  the  junction  is  abruptly  reversed  the 
stored  carriers  reach  equilibrium  after  some  storage  time  T~.  The  associ- 
ated current  may  be  defined  as 

1  "  TDg^  (3.12) 

Substituting  (3.10)  into  (3.12)  yields 

1  -  TDIsae™  -  Cs(v)^  (3.13) 

where  Cs(v)  =  I^ae01   is  the  charge-storage  capacitance.  The  charge- 
storage  capacitance  may  be  alternately  written  in  terms  of  the  dynamic 
resistance  as  C<.(v)  =  TD/r  .(v)  • 

d.  Diode  Model 

Combining  the  effects  outlined  above  results  in  the  diode 
model   in  Fig.   3.1.     The  circuit  equation  becomes 

cdM  HI   +    1d    +    4     ■      °  <3-14> 

on 


where  C.(vj  =  CD  +  C<iv). 
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Equation  (3.14)  is  in  the  form  of  (3.3)  where 


\ 


**<*»  -   -zjfr  -   yvj^  <3-15a> 


B(x(t))  =  0  (3.15b) 

The  diode  example  ( 1 1 1  o  D. )  does  not  include  parameter  temper- 
ature dependence.  The  range  of  values  studied  includes  a  spread  in 
parameter  values  as  may  be  due  to  manufacturing  tolerances,  environmental 
conditions,  and  aging.  The  Zener  breakdown  effect  and  voltage  dependence 
of  Cn  are  neglected.  Typical  parameter  values  are  used. 
2.  Modeling  for  a  p-n  Junction  Transistor 

Many  models  for  a  p-n  junction  transistor  have  been  proposed  in 
the  last  fifteen  years,  the  most  common  being  the  Ebers-Moll  model  [8,14], 
the  lumped  model  [16],  and  the  charge-control  model  [18,19].  It  has  been 
shown  that  each  model  yields  the  same  result  in  the  solution  of  a  large- 
signal  transient  problem  [20].  The  Ebers-Moll  model  is  used  in  this  study 
and  is  summarized  in  this  section. 

a.  Current-Voltage  Relationship 

The  current-voltage  equations  for  the  Ebers-Moll  model  are 
given  in  (3.16)  and  (3.17). 

=  jVcO (e^C.T)  .     JE0    (eaVE-l)       (3.16) 

E     I  -  ouaM  I  -  ajOu. 

,   =  !«!i0_  (e-'E.D  .  ^£^_  (e«vC-l)       (3.17) 

C     I  -  cuou,  I  -  ajOu, 


where 


a   =  constant  having  a  value  between  1  and  2 
aN  =  common-base  current  in  the  normal  mode 
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ou   =   common-base  current  gain  in  the  inverse  mode 

Ip0   =   saturation  current  of  collector  junction  with  zero 
emitter  current 

Ipfi   =   saturation  current  of  emitter  junction  with  zero 
collector  current 

b.  Transistor  Model 

Figs.  3.2a  and  3.3a  are  representations  of  (3.16)  and  (3.17) 
for  a  n-p-n  and  p-n-p  transistor  respectively.  Figs.  3.2b  and  3.3b  repre- 
sent the  corresponding  equivalent  circuits  if  each  diode  is  modeled  as 
outlined  in  the  preceding  section. 

The  reverse  saturation  currents  are  as  follows: 

Ir-  =  ■  ai  C0  (3.18a) 

CS     1  -  cuaN  v     ' 

IES  =  ^ECL  (3.18b) 

1  -  ajaN 

If  the  positive  senses  of  the  currents  are  as  indicated  in 

Figs.  3.2  and  3.3,  IFq  and  Ip~  are  positive  quantities  for  n-p-n  transistors 

/  x       didc 
and  negative  quantities  for  p-n-p  transistors  and  C~r(v)  =  1^      ,. 


H  i 

and  Cc;F(v)  =  TnF    de  .   TRF  is  approximated  by  the  transistor  rise 

dt 

time  and  Tnr  is  approximated  by  the  sum  of  the  transistor  storage  time  and 
the  transistor  fall  time.  The  rise,  fall  and  storage  times  are  complex 
functions  of  the  base  current,  collector  current,  ou,,  ou,  and  the  transistor 
current-gain  cut-off  frequency,  and  are  defined  in  Fig.  3.4  [17]. 

If  Cdc(v)  =  CDC  +  CSC(V)  and  Cde(v)  =  CDE  +  CSE(v)  the 
equations  for  Fig.  3.2b  become 

dvr  vr 

Cdc<v>dT-      +     '*    +     R^      "      Vde      ■      °  <3J9a' 
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Fig.  3.2.  Model  and  equivalent  circuit  for  the   n-p-n 
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Fig.  3.3.  Model  and  equivalent  circuit  for  the  p-n-p 
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Fig.  3.4.  Transistor  switching  times 
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dvF  vF 

Cde(v>  WT    +      'de    +    JT 


SHE 


which  are  in  the  form  of  (3.3)  where 


f(x(t))  = 


awi 


de 


dc 


^^    '     ^^ 


de 


ardc 


V'dc   =   ° 


1 
vc  ; 

RSHCCdc(v) 


RSHECde( 


7j) 


(3.19b) 


(3.20a) 


B 


(3.20b) 


Similar  equations  can  be  written  for  Fig.  (3.3b) 

The  transistor  example  (III.  E.)  does  not  include  the  effects 
of  parameter  temperature  dependence  and  Cp.r,  Cp.~,  TnF  and  T~-  are  assumed 
to  be  constant.  Typical  parameter  values  are  used. 


D.  DIODE  EXAMPLE 

A  diode  problem  of  considerable  interest  is  now  solved  in  order  to 
illustrate  the  modeling  and  the  mathematical  formulation  discussed  in  this 
chapter.  The  problem  considers  the  effect  of  a  radiation  current  pulse 
injected  across  a  p-n  junction  as  indicated  in  Fig.  3.5,  where  the  radiation 
current  pulse  is  defined  by 


di  (t) 

i   (t)  +  TD  -M =  A  Y(t) 

dpx  ;     R   dt  M  ' 


(3.21) 


where  y(t)   is  the  unit  pulse  given  in  Fig.   3.6. 

Replacing  the  diode  in  Fig.   3.5  by  the  model   described  above  yields 
Fig.   3.7.     Typical   values  are 

CD  =  1.0  pF 
a  =  1/.026  V" 
A  =   10.0  mA 


I     -   1.0  pA 


,-1 


R$H  =   100.0  M^ 


TD  =   l.Oys 
TR  =   0.2ys 
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r(t) 


Fig.  3.  5.  Diode  circuit. 


time  ifj.s) 

Fig.  3.6.  Radiation  current 
pulse. 


Fig.  3.7.      Equivalent  circuit 
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1 .  Analytical  Solution  for  Radiation  Current 

Equation  (3.21)  can  be  solved  using  Laplace  transforms. 


I     (s) 
PP 


A'(  -f-    "     s  +]1/T     }  °  <  t  <  10"7 

s  s        !/ J  R  (3.22) 

I    -394A.(rl  )  t>10"7 

K 


The  resulting  time  response  is  given  in  (3.23) 


0  <  t  <  10"7 


f  A.(l  -  e"t/TR) 
ipp(t)  -i  (3.23) 

I  .394  A  e-f*"10  3/TR       t  >  10"7 


2.  Piecewise-Linear  Diode  Model 

If  the  diode  current- voltage  relationship  is  assumed  to  be  piecewise 
linear  (Fig.  3.8)  i  .  =  0  for  v<  0  and  v  =  i  ,R,  for  v  >  0,  where  R,  =  Rd(v) 
=  r.(v).  The  storage  capacitance  becomes  a  constant,  C<.  =  D/R,  ,  for 
v  >  0  and  an  open  circuit  for  v  <  0.  Fig.  3.7  can  be  redrawn  as  in  Fig. 
3.9  where  RT  is  the  parallel  combination  of  R,  R«.„  and  R,  and  C  =  C<.  for 
v  >  0,  and  RT  is  the  parallel  combination  of  R  and  R<-H  and  C  =  Cp.  for  v  <  0. 

The  rise  and  fall  time  can  be  determined  using  Laplace  transforms. 

For  0  <t  <  10"7 

Z(s)  =  —^ (3.24a) 

s  +  1/RC 

Ks)  =  Ipp(s)  -^  (3.24b) 

rt  ls     s  +  1/TR;    s 
The  voltage  takes  the  form 


v(t)     =     C]   +  C2e"t/TR  +  C3e"t/RC  0  <  t  <  10"7     (3.25) 
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slope  = 


Fig.  3.8.  Piecewise  linear  model. 


Fig.  3.9.  Diode  circuit  representation. 
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For  t  >  10"7 

IS     =      .394  A.(s.\/Tr)     ■ 

E/R 
s 

(3.26) 

The  resulting  voltage  response  is  given  by 

v(t)     =     C4  +  C5e"t/TR  +  C6e"t/RC  t  >  10"7  (3.27) 

The  differential   equation  becomes 

c  3T    "    V  "  Gv  "  E/R  (3'28> 

where     G  =  1/R,. 

Equation  (3.28)  is  solved  using  trapezoidal  integration.  Several 
values  of  R,  are  used  and  the  output  voltage,  v  =  -v,  is  indicated  in 
Fig.  3.10.  The  computer  program  is  presented  in  Appendix  B.  The  solution 
in  Fig.  3.10  is  of  little  practical  value  because  of  the  wide  variation  of 
the  voltage  response  with  assumed  values  of  R,. 
3.  Taylor  Series  Approximation 

Equation  (3.28)  may  be  rewritten  as 

C(v)  dv  =  i   dt  -  G(v).v  dt  -  E/R  dt  (3.29) 

where  C(v)  =  CD  +  C$(v) 

G(v)  =  *  +  i+^r 

Using  trapezoidal  integration  for  (3.29)  yields 

[C(v(n))  +  C(v(n-1))]  ^    =     ^  [ipp(n)  +  ipp(n-l )  -  ^  ] 

(3.30) 
-  ^"  tG(v(n))-v(n)  +  G(v(n-1 ))-v(n-l )] 
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10.  ■ 
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.3 

-Rd  =  .8k 

1.2 

.6 

.9 

"*-Rd=2k 

* 

Rd  =  4k 

TIME       (/ts) 

FIG.  3.10.   NONLINEAR   EXAMPLE  —  PIECEWISE  LINEAR 
MODEL   WITH    Rd  ADJUSTED  FROM 
0.1  kft    TO    4.0  kfl. 
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Expanding  C(v(n))   and  G(v(n))-v(n)   in  a  Taylor  series  and  including 
only  the  first-order  terms  yields 

C(v(n))     =     C(v(n-1))  +  |||  •  Av         =       C(v(n-1))  (3.31a) 

3V  v(n-l) 

G(v(n)).v(n)     =     G(v(n-1  ))•  v(n-l )     +     G(v(n-l))<Av     +    |^-|  -  Av-v(n-l) 

vfn-1 ) 

=     G(v(n-l)).v(n)  (3.31b) 

where  it  is  assumed    ^-1        =0      and    7—  I        =     0   . 

v(n-l)  v(n-l) 

Substituting  (3.31)   into   (3.30)   results  in 
C(v(n-l)).Av      -     *  Ci     (.}  +  1w(R.n  -  ^ ] 


^[G(v(n-l)).(v(n)   -  v(n-l))] 


(3.32) 


Solving  for  v(n)  yields 

v(n)     =     2C(v(n-l))   -  AT-G(v(n-1))  v(rM) 
2C(v(n-l))  +  AT-G(v(n-1)) 

+     AT   .    V")  +  W"-D   -_3.  0.33) 

2C(v(n-l))  +  AT  '   G(v(n-1)) 


It  is  of  interest  to  note  that  the  same  result  is  obtained  if 
(3.28)  is  integrated  over  the  interval  (n-l)T  to  nT  as  a  linear,  constant- 
coefficient  equation  with  C(v)  and  G(v)  having  the  values  C(v(n-1))  and 
G(v(n-1))  respectively. 

This  method  is  implemented  with  an  adjustable  integration  step 
size.  The  step  size  is  adjusted  to  insure  that  C(v)  and  G(v)  do  not  vary 
by  more  than  a  fixed  percentage  from  (n-l)T  to  nT.  The  resulting  output 
voltage,  v  ,  is  indicated  in  Fig.  3.11,  where  the  maximum  that  C(v)  and 
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FIG.  3.11.  NONLINEAR    EXAMPLE --DIODE    RESISTANCE  AND 

CAPACITANCE   ARE  HELD  CONSTANT   BETWEEN 

CALCULATION  POINTS.  PERCENTAGE    VARIATION 
INDICATES  THE   MAXIMUM   ALLOWED   VARIATION 

IN    C(v)  AND  G(v)   BETWEEN    CALCULATION  POINTS. 
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G(v)   are  allowed  to  vary  is  10%,  1%  and  0.1%.     The  computer  program  is 
presented  in  Appendix  B.     This  method  becomes  more  accurate  as  the  percent- 
age variation  allowed  for  C(v)  and  G(v)   is  reduced. 

4.     Diode  Solution  using  General   Nonlinear  State-Variable  Format 
The  circuit  in  Fig.   3.7  contains  a  single  state  variable,  v(t). 
The  state  equation  for  the  circuit  can  be  represented  in  the  form  of  (3.3) 

where  /.% 

-  G.v(t)     -  I     (eav^  .  !) 

f(v(t))     =  _  S        (3.34a) 

CD  +  TDI$ae 

^^     =     -     I  T  /     av(t)  <3-34b> 

CD  +  TDIsae 

u(t)     =     ipp(t)   -  E/R  (3.34c) 

where  G  =     -=r  + 


R       RSH 

For  the  single  state-variable  problem  considered  in  this  example 
(3.9)  reduces  to 

AT  [f  +  j  B   *    (u(n)  +  u(n-l))] 

v(u)  =  v(n-])  +  ,   ATr9f  ;  3(b  ;  mr, —  (3-35) 

1  "    2  hv  av"         J 

af 
where  the  partial  derivatives  may  be  obtained  from  (3.34)  and  f,  B,  -r—  , 

oV 

and  ^B  '  u(n))  are  evaluated  at  v(n-l). 

oV 

The  computer  program  for  this  method  appears  in  Appendix  B  and  the 
results  are  indicated  in  Fig.   3.12.     The  program  adjusts  AT  automatically 
so  that     Av     never  exceeds  0.01   volts.     These  results  have  been  confirmed 
by  direct  integration  of  the  nonlinear  state  equation  using  a  fourth-order 
Runge-Kutta  integration  scheme. 
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FIG.  3.12.  NONLINEAR    EXAMPLE  —  SOLUTION   USING 
TAYLOR    SERIES   METHOD. 
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5.  Presentation  of  Data 

The  results  of  several  runs  are  presented  in  Figs.  3.13  through 
3.17.  The  recovery  time,  defined  as  the  time  for  the  voltage  to  recover 
to  90%  of  its  original  value,  is  plotted  as  the  ordinate.  The  abscissa 
represents  the  parameter  value. 

In  Fig.  3.13  the  amplitude  of  the  radiation  pulse,  A,  is  varied 
from  0.265  mA  to  20.0  mA.  When  A  <_  0.265mA  the  radiation  current  pulse 
does  not  cause  the  voltage  to  drop  below  the  90%  value. 

In  Figs.  3.14  through  3.17  CD,  R,  TR,  and  Tn  are  adjusted  with 
A  =  5.0  mA  and  A  =  10.0  mA  respectively  while  the  other  circuit  parameters 
are  held  constant  at  their  nominal  values. 

E.  TRANSISTOR  EXAMPLE 

The  effect  of  a  radiation  current  pulse  injected  across  the  base- 
emitter  p-n  junction  of  a  transistor  operating  in  the  active  region  (Fig. 
3.18)  is  considered  in  this  section.  The  radiation  current  pulse  is 
defined  by  (3.21)  where  y(t)  is  the  unit  pulse  given  in  Fig.  3.19. 

1 .     Transistor  Solution  using  General   Nonlinear  State-Variable  Format 
Modeling  the  transistor  in  Fig.   3.18  by  the  equivalent  circuit 
presented  in  Fig.   3.2  results  in  Fig.   3.20  where 


av 


C 


CSC(V)  =  TDC  TCS  ae  (3-36a) 

avr 

CSE(v)  =  TDE  I£S  ae  E  (3.36b) 

ide  =  IES  (eaVE  -  1)  (3.36c) 

idc  =     ICS(eaVC-D  (3-36d) 

Typical  parameter  values  used  are 

aN  =  1.0                                        Isc  =  O.OlnA  TQE  =  l.Onsec 
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Fig.  3.14.  Diode  recovery  time  as  a  function  of   CD 
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Fig.  3.16.  Diode  recovery  time  as  a  function  of  TR 
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Fig.  3.17.  Diode  recovery  time  as  a  function  of  TD 
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Fig.  3.18.  Transistor  circuit 
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Fig.  3.19.   Radiation  current  pulse 
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Fig.  3.20.   Equivalent  transistor  circuit 
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a  j   =  0.99 
TR  =  0.4  nsec 
R  HC  =  100  m 

RSHE  =  10°  m 


SE 


'DE 


'DC 


O.lpA 
l.OpF 
l.OpF 
0.1mA 


TDC  =  O.lys 
a       =  38.46  V 


-1 


The  circuit  equations  are  obtained  by  summing  the  currents  at  the 
collector  and  emitter  nodes  yielding 


dvE        v 
Cde(v)  dTT+R^+  V-Vdc"  1pp-  1R    =     ° 

dv  vp  vf 

Cdc(v)  dr+R^+  'dc-V'de^R*^    =     ° 


where 


*r  =  r(vcc  +  vc-vE) 


C.Jv)     =    Cqr(v)  *  C 


'dc 


SC 


'DC 


Cde(v)     -    CSE(v)+CDE 


(3.37a) 
(3.37b) 

(3.37c) 


The  biasing  resistors,  R.    and  R   ,  are  calculated  for  a  desired 

biasing  voltage  V     ,  base=col lector  voltage  vr,  and  base-emitter  voltage 

dvE         dvC 
vF  under  steady-state  conditions,  by  setting        -rr—  ,     -rp-      and  i 

equal   to  zero  in  (3.37).     The  program  used  is  presented  in  Appendix  B, 

and  the  resulting  values  are  R,    =  0.35  Mfi  and  R     =  4.5  ktt  for  V      =    10.0  V, 
3  b  c  cc 

vE  =  0.6  V  and  vr  =  -3.4  V. 

Substituting  (3.37c)  into  (3.37a)  and  (3.37b)  and  solving  for 
x     =  [vE  Vp]   yields  a  result  identical  in  form  to  (3.3)  with 


f(x(t)) 


aTE 


^e^ 


Cdc<v^ 


GTC 
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rcu   Ip^   (e     C  -  1) 


'aN   !ES   (e"VE  "   ]) 
v         Cdc<v> 


IES(eavE-  1)' 


C^IVT 


Isc(e»C 


^^ 


1) 


(3.38a) 


B(x(t))     = 


^^ 


C.   (v)   R 


C.   (v)   R 
dcv    '     c 


(3.38b) 


u(t) 


i      (t) 
PP 


cc 


(3.38c) 


where     GTF  =  -p—    +  p 

iE       Kc         KSHE 


and     G 


TC       R 


1       +4-+       ] 


Rc       RSHC 


The  iterative  solution  equation  is  of  the  form  of  (3.9)  where 
Vf  and  V(Bu_(n))  are  obtained  by  performing  the  indicated  gradient  opera- 
tion on  (3.38)  . 

The  computer  program  for  this  problem  appears  in  Appendix  B,  and 
the  results  are  indicated  in  Fig.  3.21.  The  program  adjusts  AT  automatically 
so  that  the  norm,  [e|  ,  never  exceeds  0.01  volts. 


:|       F      /(AVE)2    +    (AVC): 


(3.39) 


These  results  have  been  confirmed  by  direct  integration  of  the 
nonlinear  state  equations  using  a  fourth-order  Runge-Kutta  integration 
scheme. 
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FIG.  3.21.    TRANSISTOR   RESPONSE    TO  RADIATION 
CURRENT    PULSE. 
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2.  Presentation  of  Data 

The  results  of  several  runs  are  presented  in  Figs.  3.22  through 
3.25.  The  recovery  time  is  plotted  as  the  ordinate  and  the  abscissa 
represents  the  parameter  value. 

In  Fig,  3.22  the  amplitude  of  the  radiation  pulse,  A,  is  varied 
from  zero  to  3.0mA.  When  A  <  2.0yA  the  radiation  current  pulse  does  not 
cause  the  voltage  to  drop  below  the  90%  value. 

In  Figs.  3.23  through  3.25  TDC,  TDE,  Cnc,  CnF  and  TR  are  varied 
respectively  while  the  other  circuit  parameters  are  held  constant  at  their 
nominal  values. 
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Fig.  3.23.  Transistor  recovery  time  as  a  function  of  TDE  and  TDC 
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IV.  SENSITIVITY  ANALYSIS  OF  NONLINEAR  CIRCUITS 

A.   INTRODUCTION 

Sensitivity  analysis  of  nonlinear  circuits  and  systems  is  a  relatively 
new  field.  Recently  S,  R.  Parker  [4]  has  considered  the  sensitivity  analysis 
of  general  nonlinear  circuits  utilizing  auxiliary  coupled  networks.  The 
sensitivity  analysis  of  nonlinear  circuits  utilizing  the  concept  of  state 
space  is  covered  in  the  next  section,  followed  by  a  sensitivity  analysis  of 
the  nonlinear  diode  circuit  discussed  in  Chapter  III. 


B.  SENSITIVITY-AUGMENTED  STATE  EQUATION 

The  state  equation  format  adopted  in  Chapter  III  (repeated  below)  is 
convenient  when  considering  circuit  sensitivity. 


x(t)  =  A(x(t))-x(t)  +  c(x(t))  +  B(x(t)).u(t) 


(4.1) 


x.     3x. 
As  in  the  linear  case,  the  sensitivity  function,  S   = 


a.       3  In  a.     ■ 


may  be  obtained  by  performing  the  indicated  partial   differentiation  on 
(4.1),  yielding  the  sensitivity-augmented  state  equation. 


X 

A 

°  " 

X 

+ 

c_ 

+ 

rB 

0 

U^ 

+ 

"o 

!°  ] 

—      — 

_x 

w 

0 

I 

adJ 

xP. 

-4 

BJ 

H 

> 

BsJ 

£ 

(4.2) 
where     j  =  1,2,   ...,r  and  x§.,  AD,  BD,  Ag.,  Bs-,  uJ  are  defined  in  Chapter 


II  and 


^ 


ac 


3C, 


3  In  a.       3  In  a. 

J  J 


3  In  a. 


(4.3) 


Separating  the  set  of  first-order  partial   differential   equations  for 
each  sensitivity  parameter  in  (4.2)  yields 
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xsj  =  A  xsj  +  cs.  +  Asj  x  *  B  us.  *  Bsj  u  (4.4) 

for  j  =  1  ,2,  . ..  ,r 
Equation  (4.4)  may  be  rewritten  as 

*sj  =  isj(*>*sj)  +  B(*)MSj  +  Bsj(x»xsj)u  (4.5) 

where  j  =  1 ,2,  ...  ,r  and 

f  .(x,x  .)  =  Ax.+c.+A.x 
-SJ  --sjy      -S3       -S3  S3  - 

Applying  trapezoidal  integration  to  (4.5)  yields 
^.(n)  -  ^..(n-l)  =  |ICfsj(x(n),xsj(n))  +  fsj(x(n-l)sxsj(n-l))] 
+  ^1  [B(x(n))usj(n)  +  B(x(n-1  ))usj(n-l )]  +  f   [Bsj(x(n)  ,xsj(n) )  u(n) 
+  Bsj(x(n-l),xsj(n-l))  u(n-l)]  (4.6) 

The  functions  fs-(x_(n)  ,xs.(n))  and  B  .(x_(n)  ,x  .(n))ujn) ,  where  u[n) 
is  considered  a  constant,  can  be  expanded  in  a  Taylor  series  about  the 
point  (x_(n-l )  ,x_  .(n-1 ))  and  approximated  by  the  first  two  terms,  as  in 
Chapter  III ,  yielding 

fsj(x(n),xsj(n))  =  fsj(x(n-l),xsj(n-l)) 

+  (AXTvx)fsj(x(n-l),xsj(n-l))  +  (ax^-v^j  )fsj(x(n-l )  ,xsj(n-l ))  (4.7a) 

Bsj(x_(n),xsj(n))u_(n)  =  Bsj(x(n-1 )  ,x_sj(n-l  ))u(n) 

+  (axTv)[B  ..(x(n-l).x  ,(n-l))u(n)]  +  (AxJ1-vx.i)[B..(x(n-l )  ,x  .(n-1 ) )  u(n)] 

(4.7b) 
where 

v   =  [-*-   ...    -^-]T 

X      3X-,  3XnJ 

= [_L_   .  .  ,      9    ]T 
Vxsj     3xi  .         9xn 

SJ  SJ 
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B(x(n))»u  .(n),  where  u  .(n)  is  considered  constant,  can  be  expanded 
—    — sj         ~sj 

in  a  Taylor  series  about  the  point  x_(n-l)  and  approximated  by  the  first 
two  terms  yielding: 

B(x(n))»usj(n)  =  B(x(n-1 )).  u^n)   +  (AxTVx)[  B(x(n-1  n'-u^fn)  ] 

(4.7c) 
Substituting   (4.7)   into   (4.6)  yields 

x^n)   -  ^(n-1)     =     *  ^  +  (MTVx)^  +  (^.^ 
+    4  [B'^jCn)  +  ^(n-l))  +  (AxTVx)(Busj(n))] 


+    ^  [Bsj"(u(n)  +  u(n-l))  +  (AxTvx)(Bsju_(n))  +  (M^Vxsj)(Bsju.(n))] 

(4.8) 

where  f  .,  B  .  and  B  .u(n)  are  evaluated  at  the  point  (x(n-l)), 
-sj   sj     sj-  K     v-v   /;' 

Recall  that 


Substituting  (4.9)   into  (4.8)  and  solving  for  x.(n)  yields 

x.jtn)     r    ^(n-DH-ATtl-^ftv^^Mv^tB^n)]1)^     f1 

[fsj  +  J-  C(AXTVX)^  f  B-t^-tn)  +  ^.(n-1))  f  (AxTVx)  [Bu^n)] 

+  Bcl-(u_(n)  +  u(n-l))  +  (4xTV  )[B     u(n)]}  ]  (4.10) 

SJ  a  bj 

Equations   (3.9)  and  (4.10)  can  be  solved  directly  for  _x(n)  and  x   .(n) 

s  j 

from  the  known  values  of  _x(n-l),  ujn),  u[(n-l),  u  ^(n),  and  u.(n-l).  The 
solution  yields  the  time  response  of  the  states  and  their  corresponding 
sensitivities. 
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C.  NONLINEAR  EXAMPLE 

In  order  to  illustrate  the  technique  discussed  above,  sensitivity 
analysis  of  the  diode  circuit  in  Chapter  III  is  considered.  The  state 
equation  for  the  diode  circuit  of  Fig.  3.5  contains  a  single  state  var- 
iable. Equation  (4.10)  reduces  to 

vsj(n)  -  vs.(n-l)+AT.[fs.  +^{^ii  Av  +  B-(usj(n)  +  usj(n-l)) 

9(B-u  (n))  9(B  -u(n)) 

+  M av  +  B  .-(u(n)  +  u(n-l))  +  — ^ Av}] 

9v  SJ  9v 

AT  9fsi    9(B  ,u(n))      , 

SJ  SJ 

9f  .          3(B-u  .(n))      3(B  .u(n)) 
where  f  .,  -^-   ,  B,  B  .,  1| and  1_ are  evaluated 

at  the  point  v(n-l ) . 

The  state  equation  may  be  written  in  the  form  of  (4.1)  where 

A(v(t))  = § —  (4.12a) 

CD  +  TDIsae 

Is(eaV-D 

c(v(t))  =  --2 —  (4.12b) 

CD  +  TDIsae 

and  B(v(t))  and  u(t)  are  given  in  (3.34b)  and  (3.34c)  respectively. 

The  equations  for  the  sensitivity  functions  take  the  form  of  (4.5), 
The  sensitivity  parameters  considered  are  the  diode  parameters,  Tn  and 
Cp.,  the  radiation-pulse  parameters,  A  and  TR,  and  the  circuit  parameters, 
R  and  E.     The  terms  in  (4.5)  are  presented  in  table  4.1   for  each  of  the 
sensitivity  parameters  where 
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k3     =     ak2  k4     =     I$(eaV  -  1) 

D     =     CD  +  k2 

9i  di 

The  expressions  for    — ££ —      and  PP...      are  obtained  by  performing 

9lnA                             R                       3f,      d'(B-u     (n)) 
the  indicated  partial   differentiation  on   (3.21)  and       „  J    ,  1* , 

dV  d  V 

9("B  .u(n))     8f  .        9(B  .u(n)) 

sir- ,  i  SJ   ,  and  |J are  obtained  by  performing  the 

indicated  partial  differentiation  on  the  equations  listed  in  table  4.1. 

In  addition  to  the  sensitivity  functions  the  voltage  response, 
v  =  -v  ,  is  shown  for  (1)  nominal  parameter  values,  and  (2)  a  10%  change 
in  the  sensitivity  parameter  value.  The  response  is  presented  in  Figs, 
4.2,  4.4,  4.6,  4.8,  4.10,  and  4.12.  The  computer  program  for  generating 
the  sensitivity  functions  is  given  in  Appendix  B,  and  the  resulting 
sensitivity  functions  are  presented  in  Figs.  4.1,  4.3,  4.5,  4.7,  4.9,  and 
4.11.  The  program  adjusts  AT  automatically  so  that  |e.|  never  exceeds 
0.05  volts  and  I Av 1  never  exceeds  0.01  volts,  where 


|ej|   =  /(Av)2  +  (Avsj)2      for  j  =  1,  2,  ...  ,6      (4.13) 

These  results  have  been  confirmed  by  direct  integration  of  the  non- 
linear state  equations  using  a  fourth-order  Runge-Kutta  integration  scheme, 
It  should  be  noted  that  the  method  outlined  in  this  chapter  is  approxi- 
mately twice  as  fast  as  the  fourth-order  Runge-Kutta  integration  scheme. 

D.   INTERPRETATION  AND  APPLICATION  OF  THE  SENSITIVITY  FUNCTION 

The  definition  of  the  sensitivity  function  adopted  in  this  thesis  is 

x,      3x_.        dx. 

~    "  (3a  ./a 

infinitesimally  small  changes  in  x.  and  a..  The  sensitivity  function  may 


x.  dX.  dX. 

Sa-     =     g  TTTa       =     (aa./a   )    '     As  noted»  this  definition  applies  only  for 
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be  used  to  approximate  changes  in  a  response  due  to  fractional  changes 

in  a  parameter  as  follows 

x.   Aa. 

AXi  =  Sa  (  ~T-)  (4-14) 

J     3 

In   general  this  approximation  becomes  more  accurate  as  the  percentage 
change  in  the  parameter  value  approaches  zero.  Table  4.2  represents  a 
comparison  of  results  obtained  for  the  underdamped  response  of  the  linear 
example  in  Chapter  II.  The  output  voltages,  v  ,  v,  ,  . . .  ,  v5  ,  are 
obtained  using  resistance  values  of  R,  1.1  R  ,  ...  ,  1.5  R  ,  where  R 
is  the  nominal  value.  Approximations  for  these  voltages  are  calculated 
using  (4.14);  specifically, 

w   AR. 

va  "     \     +     SR"(TT>  <4-15> 

o 

where  v  is  the  approximate  response  after  R  has  been  changed  by  AR  , 
and  v  is  the  nominal  response.  Table  4.3  presents  a  similar  comparison 
for  the  nonlinear  diode  example.  Values  are  tabulated  for  time  intervals 
when  the  sensitivity  function  does  not  vary  rapidly*  The  results  in 
tables  4.2  and  4.3  indicate  that  the  approximation  for  the  change  in  output 
voltage  improves  as  the  amplitude  of  the  sensitivity  function  approaches  a 
constant  value.  In  this  region  correct  values  are  predicted  for  variations 
as  large  as  50%.  It  is  noted  that  (4.14)  does  not  yield  useful  results 
for  parameter  variations  as  small  as  1%  when  the  sensitivity  function  dis- 
plays sharp  peaks.  In  this  region  the  sensitivity  function  may  still  be 
used  to  determine  changes  in  circuit  responses  by  assuming  that  the  asso- 
ciated parameter  variations  cause  a  time  shift  rather  than  a  voltage 
change.  Fig.  4.13a  represents  an  idealized  switching-circuit  response 
where  the  effect  of  an  incremental  parameter  variation  is  upon  the  switch- 
ing time.  Fig.  4.13b  represents  the  corresponding  sensitivity  function. 
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Time 

Sensitivity 
Function 

8V 

vo 

ATR 
TR 

Vl 

=  .1 
va 

ATR 
TR 
V2 

=  .2 

Va 

ATR 
TR 

V5 

.5 

(ys) 

9  In  TR 
(volts) 

va 

.84 

3.0 

8.91 

8.59 

8,61 

8.27 

8,31 

7.90 

7.41 

.94 

2,0 

9.37 

9.14 

9.17 

8.91 

8.97 

8.19 

8.37 

1.13 

1.0 

9.76 

9.65 

9,66 

9,52 

9.56 

9.06 

9.26 

Sensitivity 
Function 


AE      ,      AE     0     AE     c 

— F  =   o2      — p   =   .5 


Time        9v  E  E  E 


d   In  E 
(ys)       (volts) 


v^     vi     v.     vo     v.     vc     v„ 
o  a     2     a     5     a 


.85  -10.0  9.02       10,03     10,02     11,03     11.02     14.03     14.02 

.90  -10,0  9,24       10.24     10,24     11.24     11.24     14.24     14.24 

1.20  -10.0  9.83       10.83     10.83     11.83     11.83     14.83     14.53 


Table  4.3.  Voltage  comparison  for  nonlinear  diode 
example  using  (1)  actual  parameter  changes,  (2)  approximation 

using  (4.15). 
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FIG.  4.13    (a)  IDEALIZED  CIRCUIT  RESPONSE 
AND  (b)  ASSOCIATED   SENSITIVITY 
FUNCTION. 
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The  delay  time,  T,,  due  to  the  parameter  variation  is  given  as 

Td  =  t2  -  t1  (4.16} 

Applying  (4 J 4)  to  Fig.  4.13  it  follows  that 

Aa. 
Axi  =  xss  =  h  ^"t^  (4-17) 

The  area,  A,  under  the  sensitivity  function  over  the  time  interval 


t2  -  t,  is 


t2   x 
A  =   /    S.1  dt  =  h  T.  (4.18) 


Solving  for  h  in  (4.17)  and  substituting  into  (4.18)  yields  an 
equation  for  the  delay  time  T, . 

(Aa./a.)  l2      x. 

T.     =     3—±-        /         S1     dt  (4.19) 

d  Xss  t1  aj 

Equation  (4.19)  may  be  applied  to  estimate  the  delay  time  correspond- 
ing to  a  pulse- type  sensitivity  function  by  calculating  the  area  of  the 
pulse  and  substituting  into  (4.19).  This  type  of  calculation  is  performed 
for  the  waveforms  of  Figures  4.1,  4.3,  4.5,  and  4.7  for  10%,  20%  and  50% 
parameter  variations.  The  results  are  compared  with  the  actual  time 
delays  in  table  4.4.  The  area,  A  =  /    S    dt  ,  is  calculated  assuming 

ti      aj 

that  the  amplitude  spikes  are  triangular.  Times  t-,  and  tp  are  chosen  at 

x. 
those  points  where  S  1  =  .15  h  .  The  results  indicate  time  delays  can  be 

calculated  within  10%  for  parameter  variation  as  large  as  50%.  It  is 

noted  that  calculation  of  the  area  presents  the  greatest  limitation  in  the 

accuracy  of  the  delay  time. 
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V.  CONCLUSIONS 

The  technique  of  generating  sensitivity  functions  for  linear  circuits 
using  the  concept  of  a  sensitivity  model  is  shown  to  be  an  extension  of 
the  compensation  theorem.  A  set  of  iterative  equations  is  developed  to 
yield  simultaneously  the  time  responses  of  the  circuit  states  and  of  their 
corresponding  sensitivities.  An  example  illustrates  the  procedure. 

Diode  and  transistor  modeling  for  discrete  computation  based  upon 
the  Ebers-Moll  equations  is  considered.  The  models  adopted  contain  both 
nonlinear  capacitors  and  nonlinear  resistors  for  the  junctions.  A  state- 
variable  formulation  is  developed  which  results  in  a  set  of  nonlinear 
iterative  equations  for  circuits  containing  linear  elements,  diodes,  and 
transistors.  The  use  of  these  equations  is  illustrated  by  considering  the 
effects  of  a  radiation  pulse  on  a  p-n  junction,  and  the  resulting  recovery 
times  of  the  junction  are  tabulated. 

The  foregoing  nonlinear  iterative  equations  are  extended  to  generate 
simultaneously  the  sensitivity  functions  as  well  as  the  circuit  solution. 
The  preceding  example,  which  considers  the  effects  of  a  radiation  pulse  on 
a  p-n  junction,  is  then  analyzed  from  the  point  of  view  of  a  sensitivity 
analysis.  This  analysis  indicates  that  proper  interpretation  of  the  sen- 
sitivity function  results  in  accurate  estimated  for  changes  in  circuit 
responses  due  to  parameter  variations  as  large  as  50%.  When  the  amplitude 
of  the  sensitivity  function  is  constant  the  change  in  the  circuit  response 
due  to  a  parameter  variation  is  reflected  as  a  change  in  response  ampli- 
tude. However,  in  regions  where  the  sensitivity  function  displays  impulse- 
like amplitude  peaks,  the  change  in  the  circuit  response  due  to  a  parameter 
variation  is  reflected  as  a  time  delay  or  shift  in  the  response  characteristic, 
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APPENDIX  A 
PROOF  OF  (AxTV)  f  =  (VfT)T  AX 


The  matrix  equations  Ax  Vf  may  be  written  as 


(AxTV)  f  =  [  Ax] 


ax 


+  AX.  J-     + 


v2  ax, 


+  Ax 


n  9x. 


n 


r 


Ax 


Ax 


9f, 


l  ax- 


+  AX 


+  Ax 


2  ax, 

3f, 


2  ax, 


+  Ax 


W.1 


i  axr 
af, 


+  Ax 


n  axn 


af.. 


af.. 


n    .    n 
AX,  -rcr-  +  Ax.  ■= — 


1  ax- 


'2   ax, 


af 


.  +  Ax 


n  ax 


(A.l) 


Consider  this  matrix  equation 


vf1 


~l 


3x- 


[fi 


ax. 


v 


far 


af. 


sx1      ax. 
if. 


3f. 


aXo      ax, 


afi      af' 


3xn       9Xn 


5 

ax- 


Hi 

ax, 


(A. 2) 


Hi 

ax„ 


93 


Taking  the  transpose  of  equation (A. 2) and  multiplying  by  Ax  yields 


(VfT)T  Ax 


8fl   9fl 


3x,   3x, 


9fo   3f, 


9x,   3x, 


L9X] 


3x, 


3fn 


3x. 


3f, 


3x 


3f 


3x. 


Ax. 


Ax, 


Ax 


3*i 


Ax 


1  3X- 


3f, 


AX 


1  3X- 


3f. 


1  3X- 


AX 


2  3x, 


3f. 


+   AX 


2  3X, 


3f„ 


^2-  +   ax„  — 


k2  3X, 


3fi 


+   Ax 


n  3x. 


3f, 


+   Ax 


n  3x. 


+   AX 


n  3x„ 


(A. 3) 


Equating  equation(A.3)  to  equation(A.l)  results  in  equation  (A. 4),  which 
completes  the  proof. 


(AxTV)  f  =       (VfT)T  Ax 


(A. 4) 
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APPENDIX  B 


COMPUTER  PROGRAMS 


1.   SOLUTION 
CIRCUITS 


AND  SENSITIVITY  ANALYSIS  PROGRAM  FOR  LINEAR 


C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


c 
c 
c 

c 


RFAL 

DIME 

1AMII 

2GAM( 

3GAM? 

4GAM1 

5TFMP 

6XS2( 

7TFMP 

8XP(9 

9ITIT 

TIMF 

US  M 

US(1 

EOUA 

X(NT 

XSJ( 

+  G 

WHER 


*8  ITI 
NSION 
NV(2,2 
2,1),R 
1(2,1) 
2(2,1) 
3(2,2) 
2,1)  ,X 
S3(2,  1 
1C  ),YP 
LE(12) 
=  f  .C 
ATRIX 
,1)=1A 
TIONS 
) 

NT  )       = 
AM1J*( 
E 
(  I-(T 


TLE 

AI  (2, 2), AS (2, 2), AMI (2, 2), API  (2,  2), 
),PHI(2,2),AMIT2(2,2),BS(2,1), 
Sl(2,l) ,BS2(2, 1),BS3(2,  1), 
,GAM22(2,l) ,GAM23(2, 1), GAM 11(2, 2), 
,GAM13( 2,?) ,TEMP1(2,2) ,TEMP2(2,2), 
,TFMPM2,2),US(  1,  1),XS(  2.1  )  .XSK  2,  1)  , 
SM2,1),TEMPS1(2,1),TEMPS2(2,1), 
),TEMPXS(2,  1),S  1(910 )  ,S2(91C) , S3 (910) , 
(910),AS1(2,2),AS2(2,2),AS3(2,2), 


=    2*U 

C. 

OF    THE    FORM 
PHI*X( (N-l)T) 
PHI*XSJ((N-1)T) 

X(NT )+X( (N-l)T) 


+    GAM*(U(NT)*U(  (N-l)T)  ) 

♦    GAM*(USJ(NT)-HJSJ(  (N-l)T)  ) 

+    GAM2J*(U(NT)+U( (N-l)T) ) 


/2)*A) 


(  I-(T/2)*A) 


(  I-(T 
(  I-(T 
(  I-(T 

(  I-(T 

NS       =       #    0 

NU       =       #    0 

FORMAT    AND 

THE    NUMBER 

NS=2 

NU=1 

RFAD    IN    TH 

RFAD(5,1) 

READ(5,1) 

FORMAT (8E1 

CALCULATE 

AMI  = 
T=.16E-03 
Q=T/2. 
CALL  CONST 
CALL  SUB(A 
CALL  ADD(A 
CALCULATE 

AMIINV 
CALL    RECIP 
IF(KEP-2) 
WRITF(6,4) 
F0RMAT(5HK 
GC    TO    4C 
CONTINUE 
CALCULATE 
CALL    PROD( 
CALCULATE 
CALL    CONST 
READ    IN    B 
READ(5,1) 
CALCULATE 


-1 
-1 


AMI,   <I+(T/2)*A) 

AMI  INV 


API 


PHI 
=   GAM 
=   C\M1J 
=   GAM2J 


/2)*A)   (I+(T/2)*A) 

-1 
/2)*A)   (B)*(T/2) 

-1 
/2)*A)   (ASJ*(T/2) 

-1 
/2)*A)   (BSJ)*(T/2) 
F  STATES 

F  FORCING  FUNCTIONS 
DIMENSION  STATEMENTS  MUST  BE  CHANGED  WHEN 
OF  STATES  OR  FORCING  FUNCTIONS  ARE  CHANGED 


E  I  AND  A  MATRICES 

(( AI ( I, J),J=1,NS), 1=1, NS) 

(( AS( I , J)  ,J  =  1,NS)  ,1=1, NS) 

^.2) 


(  I-( T/2)A)         API 


=   (H-(T/2)A) 


(Q, AS, NS,NS, TEMPI) 
I,TEMP1,NS,NS, AMI) 
I, TEMPI ,NS,NS,API) 
-1 
=   (I-(T/2)A) 
(.000001  ,AMI ,AMI INV,KER,NS) 
2,3,2 

ER  =  2) 


PHI  =  AMIINV^API 
AMIINV, API, NS,NS,NS, PHI  ) 

AMIT2  =  AMI  INV*(T/2) 
(Q,AMIINV,NS,NS,AMIT2) 
MATRIX 
((BS(I, J) ,J=1,NU) ,I=1,NS) 

GAM  =  AMIT2*B 
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CALL  PROD( AMIT2,BS,NS,NS,NU,GAM) 

READ  IN  BS  MATRICES  AND  CALCULATE  GAM2J  =  AMIT2*BSJ 

RE AD (5,1)  ((BS1(I,J),J=1,NU),I=1,NS) 

RE  AD (5,1)  ((BS2( I, J) , J=l ,NU) , 1=1 , NS ) 

RE  AD (5,1)  ((BS3(  I  ,  J) ,J  =  1,NU) ,I=1,NS) 

CALL  PROD( AMIT2,BS1,NS,NS,NU,GAM21) 

CALL  PROD( AMIT2,BS2,NS,NS,NU,GAM22) 

CALL  PROD( AMIT2,BS3,NS,NS,NU,GAM23) 

READ  IN  ASJ  MATRICES  AND  CALCULATE  GAM1J  =  AMIT2*ASJ 

RE  AD  (5,1)  ((  ASKI  ,  J)  ,J  =  1,NS)  ,  I=1,NS) 

READ(5,1)  ((AS2(I,J),J=1,NS),I=1,NS) 

RE AD (5,1)  ((AS3(I,J),J=1,NS),I=1,NS) 

CALL  PROD( AMIT2,AS1,NS,NS,NS,GAM11) 

CALL  PROD( AMIT2,AS2,NS,NS,NS,GAM12) 

CALL  PROD( AMIT2,AS3,NS,NS,NS,GAM13) 

WRITE(6, 12 ) 

12  FORMATS 1' ,21X ,' A  MATRIX') 
WRITE (6, 14)  ( ( AS( I ,J) .  J  =  l ,NS)  ,1=1 ,NS) 

14  FORMAT(//, 2( 10X,E10.2) ,/) 
WPITE(6,15) 

15  F0RMAT(//22X,' B  MATRIX*) 
WRITE  (6,  16)  ((BS(I,J)  ,  J=1,NI)>  ,I  =  1,NS) 

16  FGRMAT(//,1(21X,E10.2,/) ) 
WPITE(6, 13) 

13  FORMAT  (//22X,  •  PHI  MATRIX') 
WPITF(6,14)  (  (PHI ( I, J)  ,J=1  ,NS), I  =  1,NS) 
WPITE(6,24) 

24  FCRMAT(//22X, • GAM  MATRIX1) 

WPITE( 6,16)  (  (GAM(  I, J)  ,J  =  1  ,NU)  ,1=1, MS) 
WRITE(6,17 ) 

17  FORMAT(//,T27, 'SENSITIVITY  MATR I CES ' , / ,T 19 , • GAM1 1 • , T49 
1, 'GAM21' ) 

WPITF(6,19  )  (  (  (GAM1K  I  ,  J  )  , J=  1 , NS ) ,  (GAM2K  I  ,  J  )  ,  J=  1  ,NU)  ) 
1,I=1,NS) 

19  FORMAT (//,2(5X,E13.2,5X),5X,E10.2,/) 
WRITE(6,20 ) 

20  FORMAT (///,T 19 f 'GAM12' ,T49,'GAM2  2' ) 

WRITE(6,19)  ( (  (GAM12(  I,J  ),J=1,NS), (GAM22( I  ,J) , J=1,NU)  ) 
1,1=1 ,MS) 

WRITE( 6,22 ) 
22  FORMAT (/// ,T 19, 'GAM13* ,T4P, 'GAM23' ) 

WRITE (6, 19)  ( ( (GAM13( I , J ) , J=l , NS ) , (GAM23( I , J ) , J= 1 , NU ) ) 
1,I=1,NS) 

WRITE(6,31 ) 

31  FORMAT ( »1'  ,//,Tl 3, 'OUTPUT'  ,T 30 , • SE NS I T I  V  I T Y'  , T50 , • SENS 
1ITIVITY' ,T70,' SENSITIVITY' ,/,T33, 'WRT  R',T53,'WRT  L',T 
273, 'WPT  C  ,T93,'TIME' ,///) 

INITIALIZE 

RE  AD (5, 26)  ( ( X S ( I , 1 ) ,  1=  1 , NS ) , ( XS 1 ( I , 1 ) , I  =  1 ,NS ) , ( XS2 ( I , 
11)  ,1  =  1  ,NS)  ,(XS3(I  ,1),I=1,NS)  ) 

RE AD (5,26)  XP (  1 ) , YP(  1 ) , S  1  (  1 ) , S2 ( 1 )  ,  S3( 1 ) 
26  FORMAT (3E1C. 2) 

1  =  1 

Q=l. 
30  WPITE(6,32)  XS ( 1 , 1  ) , XS  1  ( 1 , 1 ) , XS2 ( 1 , 1 ) , XS3 ( 1 , 1 ) , T I  ME 

32  FORMAT(5E20.5) 

CALL  CONST(Q,XS,NS,NU,TEMPXS) 

CALL  C0NST(Q,XS1,NS,NU,TEMPS1) 

CALL  CONST(Q,XS2,NS,NU,TEMPS2) 

CALL  CONST (Q , X  S3, NS, NU, TEMPS  3 ) 

CALCULATE  STATE  VECTOR 

CALL  PROD ( PH I f TEMP XStNStNSfNUf TEMPI) 

CALL  PROD(GAM,US,NS,NU,NU,TEMP2) 

CALL  ADD(TEMP1,TEMP2,NS,NU,XS) 

CALL  ADD(XS,TEMPXS,NS,NU,TEMP4) 

CALCULATE  SENSITIVITY  VECTORS 

CALL  PR0D(PHI,TEMPS1,NS,NS,NU,TEMP1) 

CALL  PR0D(GAM11,TEMP4,NS,NS,NU,TEMP2) 

CALL  PROD(GAM2  1,US,NS,NU,NU,TEMP3) 

CALL  ADD(TEMP1,TEMP2,NS,NU,TEMP1) 

CALL  ADD(TEMP1,TEMP3,NS,NU,XS1) 

CALL  PROD( PHI , TEMP S2 , NS, NS , NU, TEMPI ) 
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CALL  PP00(GAM12,TEMP4,NS,NS,NU,TEMP2) 

CALL  PROD(GAM2  2,US,NS,NU,NU,TEMP3) 

CALL  ADD(TEMP1,TEMP2,NS,NU, TEMPI) 

CALL  ADD(TEMP1,TEMP3,NS,NU,XS2) 

CALL  PROD( PHI,TFMPS3,NS,NS,NU,TFMP1) 

CALL  PR0D(GAM13,TEMP4,NS,NS,NU,TEMP2) 

CALL  PR0D(GAM2  3,US,NS,NU,NU,TEMP3) 

CALL  A DD( TEMPI fTEMP2tNStNU, TEMPI ) 

CALL  ADD(TEMPl,TEMP3fNStNUf XS3) 

1  =  1  +  1 

TIME=TIME+T 

XPU )=TIME 

YP(I)  =  XS(1  ,1  ) 

SKI  )  =  XS1(  1,1) 

S2(I  )  =  XS2(  1, 1) 

S3(I )=XS3( 1,1) 

IF(  I  .LT.9C5)     GO    TO    30 

PLCT    THE    DATA    POINTS 

READ(5,50)     ITITLE 

50  F0RMAT(6A8) 
RFAD(5,51)  LABEL1 

51  FORMAT(AA) 

CALL    DRAW { 900, XPfYPfOtO, LABEL lf  I T I TLE , 3. E- C2 , in. , 
11,0,2,2,5,6, 1,L) 
WRITE(6,52)    L 

52  FORMAT!//, IX, 'L=» , 110) 
RFAD(5,5C)     ITITLE 

CALL    DRAW(900,XP,  S3,  1  ,  C,  L  ABEL  1 ,  I  T I  TLE  ,  3.  E- 02  , 1  0. , 
12,0,2,2,5,6, 1,L) 

WPITE(6,52)    L 

CALL    DRAW { 900 f XPf S2,2,0,LABEL1,ITITLE,3.E-02,10., 
12,0,2.2,5,6, 1,L) 

WPITE(6,52)    L 

CALL    D RAW ( 900, XP,S 1,3,0, LAB ELI, I T I TLE , 3. E- C2 , 1 0. , 
12,0,2,2,5,6,1, L) 

WPITE(6,52)    L 
40    RETURN 

END 

SUBROUTINES 

SUBROUTINE  ADD  (A,B,N,V,C) 
THIS  SUBROUTINE  ADDS  TWO  NXM  MATRICES 
DIMFNSION  A(N,M) ,3(N,M) ,C(N,M) 
DO  152  1=1, N 
DO  152  J=1,M 
152  C(  I, J)=A( I  ,J)+B(I , J) 
RETURN 
END 

SUBROUTINE  SUB  (A,B,N,M,C) 

THIS  SUBROUTINE  SUBTRACTS  TWO  NXM  MATRICES 
DIMENSION  A(N,M) ,B(N,M) ,C(N,M) 
DO  152  1  =  1, N 
DC  152  J=1,M 
152  C(  I, J)=A( I ,J)-B< I , J) 
RETURN 
FND 

SUBROUTINE  CONST ( Q, A, N, M,C ) 

THIS  SUBROUTINE  MULTIPLIES  A  NXM  MATRIX  BY  A  CONSTANT 
DIMENSION  A(N,M) ,C(N,M) 
15  DO  150  I=1,N 
DO  150  J=1,M 
150  C(  I, J)=Q*A(I ,J) 
RFTURN 
END 


97 


SUBROUTINE  PROD 
C  THIS  SUBROUTINE 
C      BY  A  MXL  MATRIX 

DIMENSION  A(Nt 

DO  1  I=1,N 

DO  1  J=1,L 
1  C(  If J)=C. 

DO  151  1=1 

DO  151  J=l 

DO  151  K=l 
151  C  (  I  ,  J  )  =C  (  I 

RETURN 

END 


(A,B, 
POST 


N,M,L,C) 
MULTIPLIES 


A  NXM  MATRIX 


MI,B(M,L)fC(NfLi 


fL 

fM 


J)+A(I ,K)*B(K,  J) 


2 

10 


11 
12 
13 

14 


15 
20 
3  0 
31 


32 

33 

35 
36 
34 

40 


41 

42 
43 


50 


SUBROU 
THIS  S 
D  I  MENS 
DO  1  I 
DO  1  J 
X(  I, J) 
DC  2  K 
X(K,K) 
DO  34 
KP=0 
1  =  0. 
DO  12 
IFIZ.G 
Z=ABS( 
KP  =  K 
CONTIN 
IFIL.G 
DO  14 
Z=A(L, 
AIL, J) 
A(KP, J 
DO  15 
Z=X(L, 
X(L, J) 
X(KP,J 
IF (ABS 
I  F  <  L  .  G 
LP1=L+ 
DO  36 
IF(A(K 
RATIO= 
DC  33 
A(K, J) 
DO  3  5 
XIK, J) 
CONTIN 
CONTIN 
DO  4  3 
II=N+1 
DO  43 

s=o. 

IFIII. 

I  IP1=I 

DO  42 

S=S+A( 

X (  II  ,  J 

KFR=1 

RETURN 

KER=2 

RETURN 

END 


TINE  RFC  IP(EP,A, X,KER,N) 

UBROUTINE  TAKES  THE  INVERSE  OF  A  NXN  MATRIX 

ION  A(N,N) ,X(N,N) 

=  1,N 

=1,N 

=C. 

=  1»N 

=  1. 

L=1,N 


K=L,N 

E.ABS( AIK,L) ) )  GO  TO  12 

A(K, L)  ) 

UE 

E.KP>  GO  TO  20 

J  =  L,N 

J) 

=  A(KP,  J) 

)  =  Z 

J=1,N 

J) 

=  X(KP,  J) 

)=Z 

(A(L,L ))  .LE.EP)  GO  TO  SO 

F.N)  GO  TO  34 

1 

K=LP1,N 

,L) .EQ.O. )  GO  TO  36 

A(K,L)/A(L,L) 

J=LP1,N 

=A(K,J)-RATIC*A(Lf J) 

J=1,N 

=X(K,J )-RATIO*X(L, J) 

UE 

UE 

1  =  1, N 

-I 

J=1,N 

GE.N)  GO  TO  43 

1  +  1 

K= I  I P1,N 

II, K )*X( K,J) 

)=(X( II, J)-S)/A(  II, II) 
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2.   DIODE  SOLUTICN  USING  PIECEWISE  LINEAR  MODEL 


3 

100 

3^0 


6 
4 


DIMENS 
REAL    L 
REAL*8 
1  PIE 

ICNT=1 
IEND=5 
READ(5 
FORMAT 
TIME=0 
RFAD(5 
FORMAT 
ALP=l. 
GCLD=R 
KKK=59 
KKKIN= 
NM1=1 
N=l 

IF(TIM 
GNEW=0 
GO  TO 

1  GNEW=P 

2  N=N+1 
CURP :<N 

HGOLD+ 
IF(V(N 
G=l./R 
CAP=TD 
GO  TO 
9  G=l./P 
CAP=CD 
5    V(N)=( 

l)/(2.* 
GOLD=G 
TIME=T 
IFtN.G 
TEMT=T 
IFITEM 
NM1  =  N 
GO  TO 
20  KKK=KK 
I  I=KKK 
XP(  I  I) 
YP(II) 
TEMT=C 
NM1  =  N 
IF(KKK 
PLOT  T 
I  F ( I  EM 
IFUCN 
CALL    D 

10,5,5, 
WRITE( 
13  FORMAT 
GO  TO 

102  IFUCN 
CALL  D 

10,5,5, 
WRITE( 

103  ICNT=I 
GO  TO 

101  CALL  D 
10,5,5, 
WRITE( 
GO  TO 

104  CALL    D 
10,5,5, 

WRITE( 
106    RETURN 
END 


ION    V(600) ,CURR(600) , YP ( 905 ) , XP ( 905 ) 
ABEL/4H  / 

ITITLE(12)/«  DIODE     SOLUTION 

CEWISE    LINEAR    MODEL       T=l.E-09  •/ 


RD,NST0P,V(1) ,XP(1),YP(1),CURR(1),TEMT 
I  10, 5E10.2) 


,3)     SI ,E,CD,RSH,R,TR,TD,RADC 

(8E10.2) 

.0 

,300)     T 

(2E1^.2 

/.026 

ADC 

Q 

KKK-1 


E.LT.l.E-07)    GO    TO    1 

.0 

2 

ADC 


)=( ( 2.*TR-T) /<2.*TR+T) > *CURR ( NM1 ) + ( T /( 2.*TR+T) )* 

GNEW  ) 

Ml) .LT.O 

+1./RSH+ 

/PD 

5 

+1./PSH 


.0)    GO    TO    9 
l./RD 


(2.*C*P-G*T)*V(NM1)+T*(CURR(N)+CURR( NM1 )-2.*E/R) 

CAP+G*T) 

NEW 

IME  +  T 

E.NSTOP) 

EMT+T 

T.GE  ..15 

4 

K+l 

-KKKIN 
=  TIME 
=-V(N) 


GO    TO    20 
E-08)    GO    TO    20 


.LT.  1493 
HE    DATA 
D.EO.l) 
T.GT.l  ) 
RAW( 900, 
1,L) 
6,13)    L 


)    GO    TO    6 

POINTS 

GO    TO    104 

GO    TO    102 

XP,YP,1 ,0,LABEL1,ITITLE,3.E-07,3.,1,0,2, 


1X,«L=« ,110) 


D)    GO    TO    101 

XP, YP, 2,0, LAB EL  1,1 T I TLE , 3. E-07 , 3 . , 1,0,2, 


103 

T.EO.IENI 

RAW( 900, 

1,L) 

6,13)    L 

CNT+1 

IOC 

RAW(9->0,XP,YP,3,0,LABEL1,ITITLE,3.E-07,3.,1,0,2, 

1,L) 

6,13)    L 

106 

RAW ( 900, XP,YP, 0,0, LAB  EL  1,1  TITLE, 3. E-07, 3., 1,0, 2, 

ltLI 

6,13)    L 
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3.   DIODE  SOLUTION  HOLDING  DIODE  CAPACITANCE  AND 
CONSTANT  BETWEEN  CALCULATION  POINTS 


RESISTANCE 


11 
12 

3 

100 


300 


20 


REAL 
DIME 

1  IT  IT 
ICNT 
IEND 
READ 
FORM 
READ 
FORM 
READ 
FORM 
TIME 
CALL 
READ 
FORM 
TEMT 
ALP= 
INIT 
RD=1 
CAP= 
V(  1) 
CURR 
GOLD 
KKK  = 
KKKI 
NM1  = 
N=l 
IF(T 
GNEW 
GO  T 
GNEW 
N=N+ 
T=2. 
1F(T 
CURR 

KGOL 
G=l. 
V(N) 

l)/(2 
RD=V 
GTTE 
CTTE 
COM 
C0N2 
C0N3 
C0N4 
IF(C 
IF(G 
T=T/ 
GO  T 
GOLO 
CAP  = 
TIME 
IF(N 
TEMT 
IF(T 
NM1  = 
GC  T 
KKK  = 
I  I=K 
XPU 
YP(  I 
NV1  = 
N=l 
TEMT 
IF(K 


*8 
NS 
LE 
=  1 
=  3 
(5 
AT 
(5 
AT 
(5 
AT 
=C 
C 
(5 
AT 
=C 
1. 
IA 
.F 
CD 

7i 

=  R 
59 
N= 
1 


ITI TLE 
ION    V( 6^0) ,CURR( 600) , XP ( 905 )  , YP < 905 )  , 
(12) 


,11)     ITITLE 

(6A8) 

,12)    LABEL1 

<A4) 

,3)     SI ,E,CD,RSH,R,TR,TO,RADC 

(8F1C.2) 

.0 

ANCEL(2) 

,3C0)    T,GAM,NSTOP 

(2E1C. 2 ,110) 

.0 

/.C26 

LIZE    RD    AND    CAP 

+50 

1C. 
)=0.C 

ADC 

0 

KKK-1 


IME.LT 

=0.C 

0    2 

=  RADC 

1 

*T 

.GE.l. 

(N)=( ( 

D+GNEW 

/R+l./ 

=( (2.* 

.*CAP+ 

(N)/(S 

ST=1./ 

ST=CD+ 

=ABS(G 

=ABS(C 

=GAM*G 

=GAM*C 

0N1.LT 

AM.EQ. 

2. 

0    9 

=  GNFW 

CTTEST 

=TIME+ 

.GF.NS 

=TEMT+ 

EMT.GE 

N 

0    4 

KKK+1 

KK-KKK 

I)=TIM 

I)=-V( 

N 


.l.E-07)    GO    TO    1 


E-09)T=l.E-09 

2.*TR-T)/(2.*TP+T) > *CURR (NM1 )+(T/l 2.*TR+T) )* 

) 

RSH+l./RD 

CAP-G*T)*V(NM1)+T*(CURR(N)+CURR(NM1  >-2.*E/R) 

G*T) 

I*(EXP( ALP*V(N> )-l. ) ) 

R+l. /RSH+l./RD 

TD*SI*ALP*EXP( ALP*V(N) ) 

-GTTEST) 

AP-CTTEST) 

TTEST 

TTEST 

.C0N3.AND.C0N2.LT.C0N4)     GO    TO    7 

0.0)    GO    TO    7 


TOP)     GO    TO    20 

T 

..16E-08)    GO    TO    20 


IN 
E 

N) 


=  0    C 

KkIlT.1493 )    GO    TO    4 
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XP(1)=0.C 
YP(1 )=10. 
PLOT    THE    DATA 
IF(  IEND.EQ.l) 
IF(ICNT.GT.l) 
CALL    DRAW(900t 
10,5,5.1,L) 
WPITE(6, 13)    L 
13    F0RMAT(//,1X, 'L*1 
GO    TO     1C3 

102  IF( ICNT.EQ.IEND) 
CALL    DRAW(900,XP, 

10,5,5,  1,L) 
WRITE(6,13)    L 

103  ICNT=ICNT+1 
GO    TO    100 

101  CALL  DRAW(900, 
10, 5.5,1, L) 
WRITE(6,13)  L 
GC  TO  106 

104  CALL  DRAW( 900, 
10,5,5, 1,L) 

WRITE(6,13)  L 
106  RETURN 
END 


POINTS 

GO  TO  104 

GO  TO  102 

XP,  YP, 1,0, LAB ELI, I T I TLE , 3. E-07, 3. , 1, 


110) 


0,2, 


GO  TO  101 

YP, 2,0, LAB ELI, IT  I TLE , 3 . E-07 , 3 . , 1,0,2, 


XP,YP, 3,0, LAB ELI, ITI TLE, 3. E-07, 3. ,1,0, 2, 


XP,YP, 0,0, LABEL  1,1  TITLE, 3. E-07, 3., 1,0, 2, 


4.   SOLUTION  AND  SENSITIVITY  ANALYSIS  PROGRAM  FOR  NONLINEAR 
DIODE  EXAMPLE 


12 


REAL 
DIME 
1XSK 
2XSM 
3S6(9 
CALL 
NSTO 
READ 
FORM 
TIME 
READ 
READ 
READ 
READ 
READ 
FORM 
TEMT 
ALP= 
RSH= 
G=l. 
U0=- 
US10 
US20 
US30 
US40 
US50 
US60 
US1N 
US2N 
US3N 
US5N 
US1T 
US2T 
US3T 
US5T 
GOLD 
KKK= 
KKKI 
NM1  = 
N=l 
IF(T 


*8  ITI 
NSION 
600  ),S 
60C)  ,S 
05),TC 
CANCE 
P=3C0. 
(5,12) 
AT(A4) 
=0.0 
(5,3)T 
(5,3) 
(5,3) 
(5,3) 
(5,3) 
AT(8E1 
=  o.C 
1./.02 
l.E+08 
/R+l./ 
E/R 
=0. 
=-E/R 
=E/R 
=C. 
=0. 
=0. 
=0. 
=-E/R 
=  E/P 
=  C. 
=0. 

=US20+ 
=US30+ 
=0. 
=  RADC 
590 

N=KKK- 
1 


TLE 

ITI TLE ( 12) ,CURR(600) ,XP(905) ,V(600) .YP<905) 

1( 905 ),XS2( 600 ),S2( 905 ) ,XS3( 600) ,S 3(905 ) , 

A(905),SCUR(60C),XS5(600),S5(9C5),XS6(600), 

UR(600) 

L(2) 

LABEL1 


US2N 
US^N 


1 


IME.LT.l.E-07)     GO    TO    1 


101 


GNEW=0.n 
GC    TO    2 

1  GNEW=RADC 

2  N=N*1 
T=2.*T 

IF(T.GE.l.E-lO)    T=1.E-10 
CGN1=SI*ALP*EXP( ALP*V(NM1)  ) 
CON2=TD*CONl 
CCN3=ALP*CCN2 

CCN4=SI*(EXP(ALP*V(NM1) )-l. ) 
DEN=CD+CON2 

DEN2=DEN**2 

0EN3=DEN**3  * 

F=-(G*V(NM1)+C0N4)/DEN 

B-l./DEN 

PF=-(G+C0Nl)/DEN+(CUN3*<G*V(NMl)-i-C0N4>  )/DEN2 

PB=-CON3/DEN2 

FSl=-( <G+C0N1)*XS1(NM1) )/DEN+( ( C0N4+G*V< NM1 ) )*<CON2 
H-CCN3*XS1(MM1)  )  )/DEN2 

FS2=-( <G+CCN1)*XS2(NM1) ) /DEN+ ( CON3*XS2 ( NM1 )*(G*V(NM1)+ 
1C0N4) I/DEN2 

FS3=-( (G+CCN1 )*XS3(NM1 )-V ( NM1 ) / R ) / DEN+( CON 3*XS3( NM 1 ) * 
1(G*V(NM1)+C0N4))/DEN2 

FS4=-( (G+C0N1)*XS4(NM1) ) /DEN*( CUN3*XS4( NM1 )* ( G*V ( NM1 >♦ 
1C0N4) ) /DEN2 

FS5=-(  <G+CONU*XS5(NMl)  )  /DEN+<  (G*  V  (  NM1J  +  CON4-  )  *  ( CD+CON3 
1*XS5(NM1) ) )/DEN2 

FS6=-( (G+C0N1)*XS6(NM1 ) ) /DEN-M  CON3*XS6 ( NM1  )* ( G*V ( NM1 ) + 
1C0N4) )/DEN2 

BS1  =  -(CUN2+C0N3*XSKNM1) )/DEN2 

BS2=-(C0N3*XS2(NM1) )/DEN2 

BS3=-(CON3*XS3<NMl) )/OEN2 

BS4=-<CON3*XSMNMl ) I/OEN2 

BS5=-(CD+CGN3*XS5(NM1) J/DEN2 

BS6=-<C0N3*XS6(NM1) )/DEN2 

PFS1=-(G+C0N1  )/DENH-(C0N3*(C0N4  +  G*V(NMn  )  )/DEN2 

PFS2=PFS1 

PFS3=PFS1 

PFS4=PFS1 

PFS5=PFS1 

PFS6=PFS1 

PXFS1=-<ALP*C0N1*XS1(NM1) )/DEN*< (G+CON1 ) *< CON2+2 . *CON3 
1*XS1(NM1) )+<G*V(NMl)+C0N4)*(C0N3  +  ALP*C0N3*XSKNMl) ) )/ 
2DEN2-(2.*CON3*(G*V(NMl ) +C0N4 ) *< C0N2+C0N3*XS1 ( MMI ) ) )/ 
3DEN3 

PXFS2=-< ALP*C0N1*XS2<NM1)  ) /DEN-M  ALP*CON3*XS2( NM1 )*<G* 
1V(NM1H-CON4H-2.*CON3*XS2<NM1H<(G«-CGN1)  )/ DEN2- ( 2 .*CON3* 
2C0N3*XS2(NM1)*(G*V(NM1)+C0N4) ) /DEN3 

PXFS3=-(ALP*C0N1*XS3(NM1)-1./R)/DENM ALP*CON3*XS3( NM1 ) 
1*<G*V(  NM1)+C0N4)+2.*C0N3*XS3(NM1)*(G+C0N1)-C0N3* 
2V(NM1)/R)/DEN2-(2.*C0N3*C0N3*(G*V(NM1)+CUN4)*XS3(NM1) ) 
3/DEN3 

PXFS4=-(ALP*C0N1*XS4(NM1) ) /DEN  + ( ALP*CON3*X S4 ( NM1 ) * ( G* 
1V(NM1)+C0NM+2.*C0N3*XSMNM1)*(G+C0N1)  )  /DE  N2-  (  2.  *CON3* 
2C0N3*XSA(NM1)*(G*V(NM1 )+C0N4) )/DEN3 

PXFS5=-< ALP*C0N1*XS5(NM1) > /DEN-M ( G+CON1 ) * ( CD+2 .*CON3* 
1XS5CNM1)  )-MG*V(NMl  )+C0N4  )  *(  ALP*CON3*XS5<  NM  1)  )  )/QEN2-( 
22.*CON3*(G*V(NMl)+CON4)*(CD+CnN3*XS5(NMll ) )/DEN3 

PXFS6=-(ALP*C0N1*XS6(NM1) ) /DEN-M  ALP*CON3*XS6  (  NM1 )*(G* 
1V(NM1)+C0N4)+2.*C0N3*XS6(NM1)*(G+C0N1  )  ) /DEN2- < 2. *CON3* 
2CCN3*XS6(NMl)*(G*V(NMl)-»-CON4)  )/OEN3 

PBS1=-C0N3/0EN2 

PBS2=PBS1 

PBS3=PBS1 

PBS4=PBS1 

PBS5=PBS1 

PBS6=PBS1 

PXBSl=-(CON3*< 1.+ALP*XS1(NM1) M/DEN2-M 2. *CON 3* ( CON2+ 
lCaN3*XSKNMl) ) )/DEN3 

PXBS2  =  MALP*C0N3*XS2(NM1) ) /DEN2-M  2.*CON3*C0N3*XS2 ( NM1 ) 
1I/DEN3 

PXBS3=-(ALP*CON3*XS3(NMl) ) /DEN2M  2  .  *CON3*CON3*XS3  (  NM1 ) 
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D/DEN3 

PXBS4=-(ALP*C0N3*XS4(NM1)  )  /DEN2M  2  .  *C0N3*C0N  3*XS4(  NM1 ) 
D/DEN3 

PXBS5=-(ALP*C0N3*XS5(NM1) ) /DEN2+ ( 2 . *C0N3* ( CD+C0N3* 
1XS5(NM1) ) ) /DEN3 

PXBS6=-(ALP*C0N3*XS6(NM1)  )  /DEN2*  (  2  .  *CGN3*C0fN3*XS6  (  NM1  ) 
D/DEN3 
9    CURR(N)=( (2.*TR-T)/(2.*TR+T) ) *CURR (NM 1 ) * ( T/( 2. *TR+T ) )* 
KGOLD  +  GNEW) 

SCUR(N)=(  <2.*TR-T  )/(2.*TR+T>  )  *SCUR  (  NM  1  )  +  (  T  /(  2  .  *TR  +  T)  )* 
KGOLD+GNEW) 

TCUR(N)=( ( 2.*TR-T)*TCUR(NM1)+T*(CURR(N)*CURR(NM1) 
l-GOLD-GNEW))/<2.*TR+T) 

UN=CURR(N)-E/R 

US4N=SCUR(N) 

US6N=TCUR(N) 

UT=UN+UO 

US4T=US4N+US40 

US6T=US60+US6N 

V(N)=V(NM1  )+(T*(F+.5*B*UT)  )/<  1  .-T* ( PF  +  PB*UN) /2 .  ) 

DELV=V(N)-V(NM1) 

VNCRM=    ABS(DELV) 

IF(VNORM.GT.VMIN)     GO    TO    5 

XS1(N)=XS1(NM1 )+<T*<FSl+.5*<QELV*<PXFSH-PB*USlN+PXBSl* 
lUN)  +  B*USlT-»-3Sl*UT)  )  )  /  (  1 .  -  (  T*  (  PF  SI* PBS 1  *UN )  )/2.  ) 

DELS1=XS1(N)-XS1(NM1) 

S1N0RM=    S0PT(DELV**2+DELS1**2) 

IF(SINGRM.GT.SIMIN)    GO    TO    5 

XS2(N)=XS2(NMl)+(  T*( FS 2*. 5*( DELV*( PXF S2+PB*US2N+PXBS2* 
1UN)4-B*US2T*BS2*UT)  )  I  /  ( 1 .  -  ( T*  {  PFS2+PBS  2*UN  )  )/2.  ) 

DtLS2=XS2(N)-XS2(NMl) 

S2NORM=    SQRT(QELV**2+DELS2**2) 

IF(S2NinRM.GT.S2MIN)    GO    TO    5 

XS3(N»=XS3(NMl)+( T*< FS3+. 5*( DEL V* ( PXF S3+PB*US3N+PXBS3* 
1UN)*-B*US3T+BS3*UT)  )  )/  (  1 .-  <  T*(  PFS3  +  PBS3*UN  )  1/2.  ) 

DELS3=XS3(N)-XS3(NM1) 

S3NORM=    SQRT<DELV**2+DELS3**2) 

IF(S3NORM.GT.S3MIN)    GO    TO    5 

XS4(N)=XS4(NM1 )  +  < T*( F S4+. 5*( DELV* ( PXFS4+PB*US4N+PXBS4* 
1UN)+B*US4T+BS4*UT) ) )/( 1 .-( T* ( PFS4+PBS4*UN) )/2. ) 

DELS4=XS4<N)-XS4(NM1) 

S4N0RM=    SQRT<DELV**2  +  n>ELS4**2) 

IF(S4N0RM.GT.S4MIN)    GO    TO    5 

XS5<N)=XS5(NM1)  +  <T*(FS5+.5*(DELV*<PXFS5  +  PB*US5N«-PXBS5* 
IUN)+B*US5T+BS5*UT| ) ) / ( 1 .- ( T*( PFS5+PBS5*UN) )/2.) 

DELS5=XS5(N)-XS5(NM1) 

S5NORM=    SQRT(DELV**2+DELS5**2) 

IF(S5NORM.GT.S5MIN)    GO    TO    5 

XS6(N)=XS6(NM1  )♦{  T*  (  FS6  +  .  5*  (  DELV*  (  PXFS6-»-PB*US6N+PXBS6* 
1UNKB*US6T+BS6*UT)  )  )/(  l.-(T*<  PF  S6*PBS6*UN)  1/2.  ) 

DELS6=XS6<N)-XS6(NM1) 

S6N0RM=    SQRT<DELV**2+DELS6**2) 

IF(S6N0RM.GT.S6MIN)     GO    TO    5 

GO    TO    7 
5    T=T/2. 

GO    TO    9 
7    GOLD=GNEW 

TIME=TIME+T 

UO=UN 

US40=US4N 

US60=US6IM 

IF(N.GE.NSTOP)     GO    TO    20 

TEMT=TEMT+T 

IF(TEMT.GE  ..15E-08)  GO  TO  20 

NM1  =  N 

GO    TO    4 
20    KKK=KKK+1 

I I=KKK-KKKIN 

XP(I I )=TIME 

YP(II)=-V( N) 

SKII  )=XS1  (N) 

S2(II)=X52(N) 
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S3U  I)=XS3(N) 

SMII)=XSMN) 

S5(I  I)=XS5(N) 

S6(II)=XS6(N) 

WRITE(6,14)  TIME,CURR(N),V(N) , XS1 ( N) , XS2 ( N  I , XS3(N), 
1XS4(N) tXS5(N),XS6(N) , T,N 
14  FGRMAT(10E12.4,I5) 

NM1  =  N 

N=l 

TEMT  =  C0 

IF(KKK.LT. 1493  )  GO  TO  4 
C      PLOT  THE  DATA  POINTS 

REAO(5,ll)  ITITLE 
11  FCRMAT(6A8) 

CALL  DRAW( 900, XP , YP, 0 , 0, LABEL  1 , I TI TLE , 3. E-07 , 3 . , 
lltO,2,2,5,5,l,L) 

WRITE(6,13)  L 
13  FORMAT(//, IX, • L= • ,110) 
C      PLOT  THE  SENSITIVITY  FUNCTIONS 

READ(5,11)  ITITLE 

CALL  DRAW (900, XP , S 1 , 0 , 0 , L ABEL  1 , I T 1 TLE , 3. E- 07 , 20. , 
10,0,2,2,5,5, 1,L) 

WRITE(6,13)  L 

READ(5,11)  ITITLE 

CALL  DRAW ( 900, XP,S2, 0,0, LAB ELI, IT  I TLE , 3. E-C7, 100 . , 
15,0,2,2, 5, 5, 1,L) 

WPITEC6.13)  L 

REA0(5,11)  ITITLE 

CALL  DRAW(9DC,XP,S3,0,0,LABEL1,ITITLE,3.E-C7,100. , 
10,0,2,2,5, 5, 1, LI 

WPITE(6,13)  L 

READ(5,11)  ITITLE 

CALL  DRAW(900,XP,S4,0,0,LABEL1,  IT  I TLE, 3. E- 07, 1 00 . , 
10,0,2,2,5,5,1,L) 

WPITE(6,13)  L 

READ(5,11)  ITITLE 

CALL  DRAW (900, XP , S5, C , C, LABEL  1 , I T I TLE , 3. E-0  7 , 2. , 
13,0,2,2, 5, 5, 1,L) 

WRITE(6,13)  L 

READ(5,11)  ITITLE 

CALL  DRAW ( 90C , XP , S6 ,0 , 0 , L ABEL  1, IT ITL E , 3. E-07, 1 0. , 
11,0,2,2,5,5,1,L> 

WPITE(6,13)  L 
106  RETURN 

END 

5.   TRANSISTOR  BIASING  PROGRAM 

CALL  CANCEL12) 
C      READ  IN  THE  TRANSISTOR  PARAMETERS 

RE  AD (5, 2)     RSHE,CDE,TDE,SIE,ALPE,ALPI 
RE AD (5, 2)     RSHC  ,CDC  ,TDC ,S IC, ALPC, ALPN 
C  READ    IN    THE    DESIRED    BASE-EMITTER,     BASE-COLLECTOR    AND 

C  BIASING    VOLTAGES 

READ(5.2)     E,X1,X2 
2    FORMAT(6E10.2) 

GC=(X1/PSHE-ALPI*SIC*(EXP(ALPC*X2)-1. >+SIE*( EXP(ALPE* 
1XD-1.  >)/(E+X2-Xl) 
RC=1./GC 

GB=(GC*(X1-E-X2)-X2/RSHC+ALPN*SIE*(EXP(ALPE*X1)-1.) 
1-S IC* ( EXP ( ALPC *X2)-1.) )/X2 
RB=1./GB 
C      WRITE  THE  VALUES  OF  THE  BIASING  RESISTORS 
WRITE(6,5)  RB,RC 
5  F0RMATC1'  ,//2E15.5) 
RETURN 
END 
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6.        SOLUTION    PROGRAM    FOR    NONLINEAR    TRANSISTOR    EXAMPLE 

REAL*8    ITITLE 

DIMENSION    F<2,1),QER<2,2),AI(2,2),TEMP(2T2). 
1V(9C5)  ,XPC9351 ,XS(2,1I  ,  TEMPI  (2  ,  2  )  ,  DERI  NV  (  2  ,  2)  , 
2B<  2,2) ,U< 2,1) ,BDER( 2,2), UT< 2,1)  ,TEMP2( 2,2) ,UTEMP(2,1) 
3, ITITLE( 12 ) 

CALL    CANCELC2) 

NS  =  2 

NU  =  1 

RE  AD  (5,20)     (<  Aid  ,  J)  ,  J  =  1,NS)  ,I  =  1,NS) 
20    FORMAT(8E10.2) 

PEAD(5,2)  RSHE,CDE,RB  ,  S  I  E, ALPE, AL P I 

READ(5,2)  RSHC,CDC,RC  ,SI C, ALPC , ALPN 

READ(5,2)  E,TDC,TDE,T,TR,RADC 
2  FORMAT(6E10.2) 

M=l 

TI=0. 

XMIN=.01 

CMIN=l.E-06 

TMAX=l.E-08 

HCNT=1. 

GE=1»/RC+1./RSHE 

GC=1./RC+1./RSHC+1./RB 

Xl  =  .6 

X2=-3.4 

V(  1)=X1-X2 

XP(1)=TI 

XSU,1)  =  X1 

XS(2,1)=X2 

UTEMP(  1,11=0. 

UTEMP(2,1)=E 
6  IFITI.GE.l .5E-06)  RADC=0. 

T=2.*T 

IF(T.GT.TMAX)     T=TMAX 
4    CCNE1=SIE*ALPE*EXP(ALPE*X1 ) 

CONC 1= SIC* ALPC *EXP(ALPC*X2) 

CONE2=TDE*CONEl 

CONC2=TDC*CONCl 

CONE3=ALPE*CONE2 

CCNC3=ALPC*CONC2 

CGNE4=SIE* (EXP(ALPE*X1)-1. ) 

C0NC4=SIC*(EXP(ALPC*X2)-1. ) 

DENE=CDE+CCNE2 

DENC=CDC+CCNC2 

F(l,  1)=<-GE*X1+X2/RC+ALPI*C0NC4-C0NE4)/DENE 

F(2,1)=(-GC*X2+X1/RC+ALPN*CGNE4-C0NC4)/DENC 

B(  1, 11  =  1. /DENE 

B(  1,2)=1./(QENE*RC) 

B(2,1)=C. 

B( 2,2)=-l./<DENC*RC) 

U(2,1)=E 

DER( 1, 1)=-<GE+C0NE1+C0NE3*F( 1,1) )/DENE 

DER(  1,2)  =  (  l./RC+ALPI*CONCl)/DENE 

DER(2, 1)=( 1./RC+ALPN*C0NE1 )/DENC 

DFR(  2,2)=-(GC+CONCH-CONC3*F(2,l  )  )/DENC 

BDER(1,21=C. 

BDER(2,1)=C. 
8  U(  1,  1)=(  (2.*TR-T)*UTEMP(l,m-2.*T*RADC)/<2  .*TR+T) 

BDER(1,1 )=-(CONE3*(U( 1, 11+E/RC) )/(DENE**2) 

BDER(2,2)=(CONC3*E)/(RC*(DENC**2) ) 

CALL  ADD  ( BDER,DER,NS,NS, TEMPI) 

G=T/2. 

CALL  CONST  ( Q , TEMPI , NS ,NS ,TEMP ) 

CALL  SUB  ( AI , TEMP, NS,NS, TEMPI) 

CALL  RECIP  ( .000001, TEMPI, DER INV, KER, NS) 

Q  =  T 

CALL  CONST   (Q,DERINV,NS,NS,TEMP) 

CALL  ADD  (U,UTEMP,NS,NU,UT) 

Q=.5 

CALL  CONST  ( Q, UT, NS,NU, UT ) 
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c 
c 
c 
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CA 
CA 
CA 
CA 
PR 


DE 
DE 
DE 
XN 
CN 
IF 
IF 
GO 
T= 
GO 
TI 
XI 
X2 

Q= 
CA 
CA 
TE 
IF 
HC 
N  = 
XP 
V( 
IF 
RE 
FO 
RE 
FO 
CA 
10, 
WR 
FO 
RF 
EN 


LL  P 
LL  A 
LL  P 
LL  A 
OCED 

VE( 

VC( 
FINE 
LX1  = 
LX2  = 
ORM  = 
ORM= 
(CNO 
(XNO 

TO 
T/2. 

TO 
=  TI  + 
=  TEM 
=  TEM 
1. 

LL  C 
LL  C 
MT=H 
(TI. 
NT  =  H 
N+l 
(N)  = 
N)=X 
(N.L 
AD(5 
RMAT 
AD(5 
RMAT 
LL  D 
0,2, 
ITE( 
RMAT 
TURN 
D 


ROD  (B,UT,NS,NS,NU,TEMP2  ) 

DD  ( TEMP2,F,NS,NU, TEMPI) 

ROD  (TEMP, TEMPI, NS,NS,NU. TEMP2 ) 

DD   <TEMP2,XS,NS,NU,TEMP ) 

URE  GOOD  FOR  SMALL 

N)-VE(N-l) 

N)-VC(N-l) 

A  NORM  =  SQRT(DELXl**2-»-DELX2**2) 
TEMP(i,l)-XS(l,l) 
TEMP(2,1)-XS(2,1) 
SQRT(0ELX1**2+DELX2**2) 
ABS(U(1, 1)-UTEMP( 1,1)) 
RM.GT.CMIN)  GO  TO  5 
RM.GT.XMIN)  GO  TO  5 
7 


8 
T 

P(  1, 
P(2, 


1) 
1) 


ONST  (Q,U,NS,NU, 
CNST  (Q,TEMP,NS, 
lCNT*l.E-08 
LT.TFMT)  GO  TO  6 
CNT+1. 


UTEMP) 
NU,XS) 


TI 

1-X2 

T.9C 

til) 

(6A8 

tl2) 

(A4) 

RAW( 

2,5, 

6,13 

(//, 


3)  GO  TO  6 

ITITLE 
) 

LABEL1 

900,XP,V  ,0,0, LAB  ELI, I T  I  TLE , 2 . E-06 , 1  .  , 

5,1,L) 

)  L 

1X,«L=» ,110) 
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